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ABSTRACT 

A  numerical  algorithm  was  developed  for  computing  the 
nonlinear  dynamic  response  of  a  ring-stiffened,  nearly  cir- 
cular cylindrical  shell  of  finite  length  under  transient, 
axisymmetric  radial  loads  of  arbitrary  axial  distribution. 
Nonlinear  Donnell-type  equations  were  solved  using  Fourier 
series  expansions  of  the  dependent  variables  in  the  circum- 
ferential coordinate,  modified  finite  difference  approxi- 
mations of  the  axial  derivatives,  and  Newmark ' s  beta-method, 
combined  with  Gauss  elimination,  for  the  time  integration. 

The  response  of  a  simply  supported  shell  under  an 
exponentially  decaying,  uniform  pressure  was  computed  for 
peak   pressures  and  total  impulses  between  the  static 
buckling  limit  and  the  dynamic  buckling  limit.   Near  the 
dynamic  buckling  limit,  the  exponential  growth  of  the  static 
buckling  modes  dominated;  but  as  the  peak  pressure  was 
reduced,  the  parametrically  excited  Mathieu  modes  became 
increasingly  important.   The  significance  of  damping,  the 
initial  imperfections,  and  nonlinear  coupling  was  also 
investigated. 
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CHAPTER  I 


INTRODUCTION 


I.   SCOPE  OF  THE  INVESTIGATION 

This  dissertation  considers  the  general  nonlinear 
dynamic  response  of  a  nearly  circular,  nearly  cylindrical, 
ring-stiffened  shell  of  finite  length  subjected  to  a  tran- 
sient, axi symmetric  radial  load. 

The  research  was  part  of  a  project  for  the  Naval 
Weapons  Laboratory,  Dahlgren,  Virginia,  to  analytically 
predict  the  dynamic  response  of  the  large  NWL  shock 
tube  [1]   under  operational  loading  conditions.   The  various 
sections  of  the  shock  tube  can  be  idealized  as  thin,  ring- 
stiffened,  circular  cylindrical  shells  having  either  clamped 
or  simply  supported  ends.   The  operating  loads  can  be  de- 
scribed as  internal,  propagating,  axisymmetric  pulses  having 
pressures  both  higher  and  lower  than  ambient  atmospheric 
pressure.   In  the  unloaded  condition,  the  midsurface  of  the 
shell  deviates  from  a  perfect  circular  cylinder,  but  the 
unloaded  shell  was  assumed  to  be  stress  free.   Due  to  this 
initial  imperfection  of  the  shell,  the  dynamic  response  will 
not  be  axisymmetric  even  though  the  loading  may  be  axi- 
symmetric . 


^Numbers  within  brackets  refer  to  items  in  the  List 
of  References,  page  179. 


23 


Although  there  has  been  a  great  deal  of  recent 

2 
interest  in  problems  of  this  nature,   there  are,  in  fact, 

no  previous  analyses  that  are  fully  capable  of  solving  the 
problem. 

Consequently,  an  algorithm  was  developed  for  the 
numerical  solution  to  the  Donnell-type ,  nonlinear  differen- 
tial equations  that  describe  the  general  dynamic  response 
of  the  idealized  shock  tube  sections  to  a  transient,  axi- 
symmetric,  radial  load  of  arbitrary  time  history  and  axial 
distribution.   Some  of  the  features  of  the  method  were 
demonstrated  by  computing  the  dynamic  response  of  a  simply 
supported  shell  under  an  exponentially  decaying,  uniform 
pressure.   The  peak  pressures  and  the  total  impulses  of  the 
loads  considered  fall  between  the  static  buckling  limit  and 
the  dynamic  buckling  limit.   In  addition,  the  significance 
of  viscous  damping,  of  the  initial  imperfections,  and  of 
the  nonlinear  coupling  was  investigated. 


2 
The  analysis  is  applicable  to  structures  other  than 

the  NWL  shock  tube.   For  instance,  a  missile  or  submarine 

engulfed  by  a  blast  front  presents  a  similar  structural 

dynamics  problem. 


24 


II.   HISTORICAL  REVIEW 

Rings  and  Unstiffened  Cylinders 

The  nonlinear  dynamic  response  of  circular  cylindrical 
shells  has  received  considerable  attention  over  the  past 
twenty  years,  with  the  Russians  doing  much  of  the  early 
work  in  the  field.   For  example,  the  1956  book  by  the 
Russian  V.  V.  Bolotin  [2]  is  a  classic  in  the  field  of 
dynamic  instability.   In  it  the  works  of  most  of  the  early 
Russian  investigators  are  described  in  some  detail. 

One  of  the  earliest  articles  to  appear  in  the  American 
literature  was  a  1961  translation  from  the  Russian  of  a  1958 
paper  by  V.  L.  Agamirov  and  A.  S.  Vol'mir  [3].   In  this 
paper  the  Bubnov-Galerkin  method  was  used  to  obtain  an 
approximate  solution  for  the  response  of  a  closed  cylindri- 
cal shell  subjected  to  a  hydrostatic  loading  that  increased 
linearly  with  time.   They  assumed  a  function  for  the  spatial 
variation  of  both  the  initial  imperfection  and  the  total 
deflection,  allowing  for  variation  in  the  axial  and  circum- 
ferential directions.   Thus,  the  problem  was  reduced  to  the 
numerical  solution  of  a  single  nonlinear  ordinary  differen- 
tial equation  in  the  time  dependent  amplitude  of  the  dis- 
placement function. 

In  1961,  J.  C.  Yao  [4]  published  one  of  the  first 
American  papers  to  appear  in  the  field.   He  investigated 
the  stability  of  an  infinite  cylinder  under  a  uniform  radial 
pressure  applied  as  a  rectangular  shaped  function  of  time. 
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The  motion  was  assumed  to  consist  of  a  uniform  radial  dis- 
placement plus  one  flexural  mode  having  two  waves  around 
the  circumference.   The  only  coupling  retained  was  the 
influence  of  the  axisymmetric  mode  on  the  flexural  mode. 
The  flexural  mode  was  shown  to  respond  in  an  oscillatory 
or  exponential  manner  depending  on  the  magnitude  of  the 
pressure . 

The  Agamirov  and  Vol'mir  paper  [3]  was  followed  in 
1962  by  another  translation  from  the  Russian  of  a  1960 
paper  by  Y.  I.  Kadashevich  and  A.  K.  Pertsev  [5].   They 
investigated  the  response  of  a  finite  cylindrical  shell 
loaded  by  uniform  radial  pressure  applied  as  a  step  func- 
tion in  time.   They  used  combinations  of  two  axisymmetric 
modes  and  one  nonaxisymmetric  mode  to  represent  the  initial 
imperfection  and  the  radial  displacement.   However,  these 
modes  did  not  satisfy  the  boundary  conditions.   Using 
Lagrange's  equations,  ordinary,  coupled,  nonlinear  differ- 
ential equations  were  obtained  which  accounted  for  the 
radial  inertia  of  all  three  oscillatory  modes.   These 
equations  were  integrated  numerically  to  obtain  the  time 
histories  of  the  modes.   The  critical  load  was  defined  as 
the  load  which  caused  stresses  in  excess  of  the  ultimate 
strength  of  the  shell  material.   This  paper  is  considered 
most  significant  because  it  retained  all  the  nonlinear 
coupling  between  the  three  modes.   Most  later  works  do  not 
retain  all  of  the  coupling  terms. 
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Another  very  significant  paper  was  published  by 
Yao  [6]  in  1963.   He  investigated  the  dynamic  stability 
of  closed  cylindrical  shells  under  several  combinations 
of  harmonic  and  static  loads.   He  assumed  the  radial 
motion  to  consist  of  a  uniform  quasi-static  displacement 

plus  a  small,  nonaxisymmetric ,  flexural  perturbation  that 

3 
satisfied  the  boundary  conditions.    As  in  his  earlier 

work  [4],  Yao  retained  only  the  influence  of  the  axi- 
symmetric  mode  on  the  flexural  mode.   The  governing  equa- 
tion for  the  motion  of  the  perturbation  was  an  ordinary 

differential  equation  with  a  harmonically  varying  coeffi- 

4 
cient,  namely,  the  classical  Mathieu  equation.    When  a 

flexural  mode  corresponded  to  a  point  in  an  unstable  region 

of  the  Mathieu  stability  chart,  the  response  of  the  shell 

was  assumed  to  be  unstable.   Yao  concluded  that  a  shell  not 

only  can  withstand  static  loads  up  to  the  classical 

buckling  load,  but  also  can  endure  an  additional  harmonic 

load  so  long  as  the  harmonic  component  does  not  render  any 

of  the  flexural  modes  unstable. 


3 
The  radial  inertia  of  the  axisymmetric  mode  was 

neglected. 

4  . 

The  solution  to  the  Mathieu  equation  is  either 

oscillatory  and  bounded,  or  oscillatory  and  unbounded, 

depending  on  the  parameters  [7] .   The  unbounded  condition 

is  known  as  parametric  resonance  or  dynamic  instability. 

The  distinction  between  a  bounded  response  and  one  that 

is  unbounded  is  determined  from  the  Mathieu  stability 

chart . 
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In  late  1963,  J.  D.  Wood  and  L.  R.  Koval  [8]  pre- 
sented the  results  of  their  investigation  of  the  stability 
of  an  axially  loaded  shell  subjected  to  a  time  dependent 
uniform  radial  pressure.   Their  assumptions  and  method 
of  solution  were  essentially  the  same  as  those  employed 
by  Yao  [4,  6].   However,  they  did  not  neglect  the  radial 
inertia  of  the  axisymmetric  mode  as  Yao  had  done. 

In  1964,  J.  N.  Goodier  and  I.  K.  Mclvor  [9]  published 
a  paper  which  considered  the  response  of  an  infinitely 
long  shell  subjected  to  nearly  axisymmetric,  radial  impulse 
of  low  intensity.   Mathieu's  equation  governed  the  early 
growth  of  the  flexural  modes.   Since  the  impulse  was  of 
low  intensity,  a  single  dominant  flexural  mode  could  be 
identified  from  the  Mathieu  stability  chart.   They  assumed 
the  radial  displacement  to  consist  of  a  uniform  radial 
mode  plus  the  single  dominant  flexural  mode.   The  fully 
coupled  motion  of  these  two  modes  was  computed  numerically. 
During  the  motion,  the  axisymmetric  mode  and  the  flexural 
mode  exchanged  energy  in  a  cyclic  manner  because  of  the 
total  coupling  between  the  two  modes.   Although  the 
flexural  mode  corresponded  to  an  unstable  point  on  the 
stability  chart,  its  motion  remained  bounded  since  only 
a  finite  amount  of  energy  was  imparted  to  the  shell  by 
the  load.   This  paper  demonstrated  that  the  motion  of  a 
flexural  mode  is  not  necessarily  unbounded  even  though 
its  early  growth  is  associated  with  a  point  in  an  unstable 
region  of  the  Mathieu  stability  chart. 
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H.  E.  Lindberg  [10],  in  a  companion  to  the  above 
paper,  investigated  the  dynamic  buckling  of  thin  shells 
under  very  intensive  impulsive  loads.   He  found  that  many 
short  wavelength  flexural  modes  were  excited  to  large 
amplitudes,  and  that  buckling  occurred  as  the  result  of 
hyperbolic  growth,  rather  than  any  oscillatory  growth. 

In  1964,  R.  S.  Roth  and  J.  M.  Klosner  [11]  investi- 
gated the  dynamic  response  of  a  cylindrical  shell  subjected 
to  an  axisymmetric ,  axial  step  load.   They  used  a  Donnell- 
type  shell  theory  and  assumed  that  the  radial  displacements 
and  the  initial  imperfection  were  comprised  of  a  linear 
combination  of  four  discrete  modes,  all  of  which  did  not 
satisfy  the  boundary  conditions.   Radial  inertial  forces 
in  all  modes  were  considered.   The  Galerkin-Ritz  method 
was  used  to  obtain  the  coupled,  ordinary,  nonlinear 
differential  equations  in  the  time  dependent  coefficients. 
The  resulting  equations  were  then  integrated  numerically. 
Prior  to  this  work,  the  Russians  Kadashevich  and  Pertsev  [5] 
were  apparently  the  only  ones  who  had  accounted  for  the 
simultaneous,  coupled  motion  of  more  than  one  flexural  mode. 

In  1965,  a  survey  article  by  R.  M.  Evan-Iwanowski  [12] 
reviewed  the  state  of  the  art  on  the  parametric  resonance 
of  structures.   The  bibliography  is  very  comprehensive  with 
numerous  references  to  the  early  Russian  work.   In  this 
article  Evan-Iwanowski  suggested  that  future  work  should  be 
devoted  to  study  of  the  interaction  between  the  dynamic  and 
parametric  responses. 


29 


M.  P.  Bienick,  T.  C.  Fan,  and  L.  M.  Lackman ' s  1966 
paper  [13]  investigated  the  dynamic  stability  of  a  cylin- 
drical shell  under  a  uniform  radial  pressure,  applied 
harmonically  or  as  a  ramp  function  of  time.   The  axisym- 
metric  deflection  state  varied  as  a  half  sine  wave  along 
the  axis,  thus  satisfying  the  simply  supported  boundary 
conditions.   One  flexural  mode  of  arbitrarily  small  magni- 
tude was  assumed  to  comprise  the   perturbed  state.   It 
also  satisfied  the  boundary  conditions.   The  effects  of 
the  radial  inertia  in  both  modes  were  included.   The 
stability  problem  was  reduced  to  Mathieu ' s  equation  for 
both  types  of  loads.   Under  the  ramp  loading,  they  found 
that  not  one,  but  several  flexural  modes  were  unstable 
if  the  shell  was  short. 

The  first  works  devoted  to  the  dynamic  stability  of 
a  cylindrical  shell  loaded  by  a  moving  axisymmetric  load 
were  by  G.  A.  Hegemeir  [14,  15]  in  1966  and  1967.   The 
shell  was  infinite  in  length,  and  the  load  moved  at  con- 
stant velocity,  leading  to  a  steady  state  condition. 
Hegemeir,  like  Yao  [6]  and  Bienick,  Fan,  and  Lackman  [13], 
investigated  the  stability  of  the  nonaxisymmetric  response 
with  respect  to  arbitrarily  small  perturbations.   However, 
the  governing  equation  was  not  the  Mathieu  equation. 

In  1966,  T.  C.  Fan  [16]  retained  the  axial  and  circum- 
ferential inertia  in  an  investigation  of  the  dynamic 
stability  of  a  simply  supported  cylinder  loaded  by  uniform 
radial  pressure  applied  as  a  ramp  function  of  time.   His 
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results  showed  that  neglecting  the  axial  and  circumferential 
inertial  forces  has  little  effect  on  the  stability  of  small 
flexural  perturbations.   This  seems  to  be  the  only  study 
of  the  significance  of  the  axial  and  circumferential  iner- 
tial forces. 

In  1966,  M.  H.  Lock  [17],  like  Goodier  and  Mclvor  [9], 
studied  the  response  of  a  ring,  or  an  infinitely  long 
cylinder,  but  under  step  loads  instead  of  impulsive  loads. 
Using  finite  deflection  theory,  he  computed  the  coupled 
motion  of  the  axisymmetric  mode  and  one  nonaxisymmetric 
mode.   The  conditions  for  instability  to  small  perturba- 
tions were  shown  to  be  necessary  for  the  existence  of  large 
displacements,  but  not  sufficient.   Significant  to  the 
present  investigation,  he  showed  that  a  flexural  mode  could 
be  excited  to  large  amplitude  in  either  of  two  ways:  (1) 
When  the  magnitude  of  the  load  exceeded  the  static  buckling 
load  of   a   mode,  that  mode  grew  in  essentially  an  expo- 
nential manner  with  time  until  limited  by  nonlinear  effects. 
(2)  A  flexural  mode  could  also  be  excited  to  large  ampli- 
tudes by  the  oscillation  of  the  axisymmetric  mode,  and  the 
growth  was  an  oscillation  of  increasing  amplitude,  as  in 
parametric  resonance. 

A  1965  report  by  H.  E.  Lindberg,  et  al_.  [18]  and  the 
subsequent  1968  article  by  D.  L.  Anderson  and  H.  E. 
Lindberg  [19]  are  of  considerable  interest  to  this  disser- 
tation.  They  presented  two  theoretical  models,  a  tangent 
modulus  model  for  plastic  behavior  and  an  elastic  model, 

31 


for  predicting  the  dynamic  loads  required  to  buckle  a 
cylindrical  shell.   Only  the  elastic  model  is  of  interest 
here.   This  model  made  use  of  a  DonneU-type  shell  theory 
which  accounted  for  initial  imperfections.   The  radial 
motion  was  assumed  to  consist  of  a  uniform  radial  displace- 
ment plus  nonaxisymmetric  flexural  modes  that  satisfied 
the  simply  supported  boundary  conditions.   The  spatial 
distribution  of  the  initial  imperfection  was  taken  to  be 
the  same  as  that  assumed  for  the  flexural  mode  shapes. 
Following  Yao  [4,  6],  Wood  and  Koval  [8],  and  Bienick, 
Fan,  and  Lackman  [13],  the  influence  of  the  axisymmetric 
mode  on  the  flexural  modes  was  the  only  coupling  retained. 
No  coupling  was  retained  which  might  extract  energy  from 
the  uniform  radial  mode,  nor  was  any  coupling  between 
flexural  modes  retained.   An  exponentially  decaying  uniform 
load  and  a  linearly  decaying  uniform  load  were   con- 
sidered.  The  maximum  displacement  of  several  flexural 
modes  was  found  by  numerical  integration,  and  buckling 
was  assumed  to  have  occurred  whenever  the  maximum  displace- 
ment of  any  mode  exceeded  one  thousand  times  the  initial 
imperfection  in  that  mode.   Only  the  lower  frequency, 
exponentially  growing  buckling  modes  were  considered.   The 
parametrically  excited  modes  were  not  treated.   The  authors 
reasoned  that  these  modes  would  dissipate  energy  locally  by 
plastic  flow  and,  hence,  would  eventually  dampen  out  before 
causing  deflections  large  enough  to  buckle  the  shell. 
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The  most  recent  work  is  the  paper  by  I.  K.  Mclvor  and 
E.  G.  Lovell  [20]  which  appeared  in  late  1968.   They 
studied  the  dynamic  response  of  finite  length  cylindrical 
shells  subjected  to  a  uniform  radial  impulse.   They  assumed 
a  set  of  modal  displacement  functions  that  simplified  some 
of  the  algebra,  but  did  not  satisfy  end  conditions  of 
usual  interest.   By  first  using  a  small  nonaxisymmetric 
perturbation  approach  that  retained  only  the  influence  of 
the  axisymmetric  mode  on  the  flexural  mode,  Mathieu's 
equation  was  obtained.   They  then  determined  from  the 
Mathieu  stability  diagram  the  four  or  five  flexural  modes 
having  the  highest  initial  growth  rate.   Using  these 
modes  and  the  axisymmetric  mode,  and  retaining  all 
coupling  between  modes,  they  employed  the  Runge-Kutta 
integration  method  to  compute  the  response  of  each  mode 
to  the  uniform  radial  impulse.   A  slight  irregularity 
in  the  impulse  was  used  to  initiate  the  response  of  each 
nonaxisymmetric  mode.   They  found  that  the  modes  under- 
went beat-like  oscillations  accompanied  by  a  continual 
exchange  of  energy  between  modes.   The  stress  levels 
attained  during  the  nonsymmetric  motion  were  much  higher 
than  those  developed  in  an  axisymmetric  dynamic  response 
to  the  same  load. 

Dynamic  Stability  of  Ring-stiffened  Shells 

The  treatment  of  discrete  ring  stiffeners  in  a  non- 
linear dynamic  analysis  would  provide  a  suitable  topic  for 
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a  dissertation  in  itself.   In  the  course  of  this  research, 
a  brief  study  of  the  problem  was  conducted,  and  a  seemingly 
suitable  method  of  analysis  was  selected.   It  is  not  the 
only  method,  and  perhaps  not  even  the  best  method.   A 
complete  review  of  the  related  literature  is  deemed  to  be 
beyond  the  scope  of  this  study.   Only  a  few  references 
are  mentioned  to  provide  a  starting  point  for  an  interested 
reader . 

In  a  1967  dissertation,  W.  K.  Dietz  [21]  investigated 
the  dynamic  response  of  eccentrically  reinforced,  circular 
cylindrical  shells  subjected  to  dynamic  axial  loads.   He 
considered  only  shells  with  closely  spaced  stiffeners 
and  treated  the  structure  as  an  equivalent  orthotropic 
shell.   The  Galerkin  method  was  used  to  obtain  a  set  of 
ordinary,  nonlinear  differential  equations  in  time  which 
were  integrated  by  employing  a  predictor-corrector  method. 

J.  M.  Klosner  and  H.  N.  Franklin  [22]  presented  a 
report  in  1968  which  dealt  briefly  with  the  dynamic  insta- 
bility of  long  reinforced  cylindrical  shells  under  axial 
loads.   They,  like  Dietz,  treated  an  equivalent  ortho- 
tropic  shell  and  used  the  Galerkin  method  to  obtain  the 
governing  differential  equations  of  motion.   The  equations 
were  integrated  by  the  Runge-Kutta  method. 

The  analysis  of  the  rings  contained  in  this  disser- 
tation is  essentially  an  adaptation  of  a  method  used  by 
D.  L.  Block  [23]  in  a  1968  report  dealing  with  the  influence 
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of  discrete  ring  stiffeners  and  prebuckling  deformations 
on  the  static  buckling  of  axially  loaded  shells. 

III.   NEW  FEATURES  OF  THE  INVESTIGATION 

The  studies  which  have  been  reviewed  in  the  proceding 
section  appear  inadequate  for  one  or  more  of  the  following 
reasons : 

1.  The  motion  can  not  be  predicted  after  onset 

of  nonaxisymmetric  motion. 

2.  The  boundary  conditions  are  not  satisfied 

by  the  assumed  displacement  functions. 

3.  Coupling  due  to  nonlinear  effects  is  only 

partially  retained. 

4 .  The  assumed  displacements  are  not  general 

enough;  or  an  insufficient  number  of  modes 
have  been  included  to  encompass  both  the 
higher  frequency,  parametrically  excited 
modes  and  the  exponentially  growing,  lower 
frequency,  buckling  modes. 

5.  The  analyses  do  not  handle  loads  which  depend 

on  the  axial  coordinate. 
The  algorithm  developed  in  this  dissertation  is  more 
general  than  any  of  the  earlier  analyses  and  does  not  con- 
tain any  of  the  above  shortcomings.   This  is  accomplished 
by  the  use  of  the  following  new  features: 
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1.  The  displacements  are  expanded  in  Fourier 

series  in  the  circumferential  coordinate. 
The  series  are  then  truncated,  and  only 
the  significant  modes  are  retained. 
Those  modes  which  do  not  respond  can  be 
discarded.   As  many  as  fifteen  modes  are 
included,  and  this  number  can  be  increased 
at  the  expense  of  increased  computing  time. 
In  this  study  all  of  the  nonlinear  coupling 
between  the  modes  is  retained,  and  enough 
modes  are  included  to  simultaneously  en- 
compass the  parametrically  excited  modes 
and  the  exponentially  growing  buckling 
modes . 

2.  A  set  of  modified  finite  differences  is 

employed  for  the  derivates  with  respect  to 
the  axial  coordinate  for  the  first  time 
in  a  nonlinear  shell  problem.   These 
modified  differences  are  shown  to  be 
superior  to  the  conventional  difference 
approximations  to  derivatives. 

3.  The  iterative  Newmark  beta-method  is  employed 

for  the  time  integration  of  the  radial 
equations  of  motion.   The  axial  and  circum- 
ferential displacements  are  obtained  by 
Gauss  elimination  during  each  iteration. 
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4.  The  boundary  conditions  are  consistently 

observed . 

5.  Axisyminetric  radial  loads  of  arbitrary  time 

variation  and  axial  form  can  be  readily 
handled . 

6.  Initial  imperfections  symmetric  about  a 

plane  through  the  axis  and  of  arbitrary 
axial  distribution  are  accounted  for. 

7.  Viscous  damping  forces  are  included  in  the 

equations  of  motion. 

8.  The  effects  of  discrete,  eccentrically 

located  ring  stiffeners  have  been  included 
in  the  analysis. 

IV.   ORGANIZATION  AND  PREVIEW 

The  main  body  of  the  dissertation  is  divided  into  six 
chapters.   The  present  introductory  chapter  is  followed  by 
the  development  of  the  equations  of  motion  in  Chapter  II. 
In  Chapter  III,  a  numerical  model  is  developed  that  replaces 
the  partial  differential  equations  of  Chapter  II  with 
difference-differential  equations  in  the  time  variable.   In 
Chapter  IV,  methods  for  generating  a  solution  to  the 
difference-differential  equations  are  considered,  the  com- 
puter program  is  briefly  described,  and  the  tests  which 
were  conducted  to  verify  the  computer  program  are  discussed. 
Numerical  results  for  a  problem,  a  problem  similar  to  one 
considered  earlier  by  Anderson  and  Lindberg  [19]  and  Mclvor 
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and  Lovell  [20],  are  presented  in  Chapter  V  along  with 
studies  of  the  significance  of  damping,  of  the  initial 
imperfection,  and  of  nonlinear  coupling.   Chapter  VI 
contains  a  summary  of  the  work,  the  new  findings,  and 
recommendations  for  future  work  in  the  field. 

Four  appendices  contain  some  of  the  minor  details 
of  the  development  of  the  algorithm  and  a  complete  docu- 
mentation of  the  computer  program,  including  instructions 
for  preparation  of  input  data  cards  and  other  information 
necessary  for  the  operation  of  the  program. 


CHAPTER  II 

FORMULATION  OF  THE  PROBLEM 

I .   INTRODUCTION 

In  this  chapter  the  governing  nonlinear  differential 
equations  for  the  dynamic  response  of  a  ring-stiffened, 
nearly  circular  cylindrical  shell  of  finite  length  are 
derived.   A  Donnell-type  nonlinear  theory  [24,  25]  that 
accounts  for  initial  imperfections  in  the  shape  of  the  shell 
is  used  to  describe  the  shell  response;  the  discrete 
circular  rings  are  treated  by  linear  beam  and  elementary 
torsion  theories.   The  shell  and  ring  materials  are 
assumed  to  be  homogeneous,  isotropic,  and  elastic;  and  the 
rings  may  be  fabricated  from  a  different  material  than  the 
shell.   The  centroids  of  the  ring  cross  sections  are  eccen- 
trically located  with  respect  to  the  midsurface  of  the 
shell.   Radial  inertial  forces  of  the  shell  and  the  rings 
have  been  included;  but  their  axial,  circumferential,  and 
rotatory  inertial  forces  have  been  neglected.   Viscous 
type  damping  has  been  included  to  account  for  dissipative 
effects . 

II.   COORDINATE  SYSTEM 

The  motions  of  the  shell  of  length  L  will  be  referred 
to  the  fixed  coordinate  system  shown  in  Figure  2.1.   The 
positive  direction  of  the  axial  coordinate  x,  the 


39 


Midsurface   of 
Perfect   Shell 


FIGURE        2.1 

GEOMETRY   AND  COORDINATES 
OF    SHELL 
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circumferential  coordinates  y  and  6 ,  the  radial  coordinate 
z,  the  axial  displacement  U,  the  circumferential  displace- 
ment V,  and  the  radial  displacement  W*  are  as  indicated  in 
Figure  2.1.   The  origin  of  the  coordinate  system  is  located 
at  the  midsurface  of  a  perfectly  circular  cylinder  of 
radius  a.   The  ends  of  the  shell  are  located  at  x  =  0  and 
x  =  L.   The  shell  thickness  is  denoted  by  h. 

The  motion  of  a  ring  will  be  described  by  the  dis- 

placements  of  the  centroid  of  the  ring  cross  section  U  , 

r       r 
V  ,  and  W  ,  and  the  rotation  about  a  normal  to  the  cross 

section  $,  with  positive  directions  as  indicated  in  Figure 

2.2.   The  distance  e  between  the  centroid  of  a  ring  cross 

section  and  the  midsurface  of  the  shell  is  positive  when 

the  rings  are  on  the  outside  of  the  shell. 

III.   STRAIN-DISPLACEMENT  RELATIONS 

Shell 

In  the  Donnell-type  nonlinear  theory,  the  relations 
between  the  midsurface  strains  and  displacements  are  given 
by  [25] 

2       _ 
e   =  U,   +  h    (W,  )   +  W,  W,  (2.1a) 

XX  X  XX 

W  2        — 

e   =  V, +  3§  (W,  )   +  W,   W,  (2.1b) 

y      y    a    2     y        y    y 


y         =    U,   +  V,   +  W,   W,   +  W,  W,   +  W,  W,       (2.1c) 
'Xy     'y     'x     'x   'y     'x'y     'x'y 
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^E 
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FIGURE  2.2 

DISPLACEMENTS   AND  ROTATION 
OF    A       TYPICAL     RING 
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Xx  =  W,xx  (2. Id) 


X   =  W,  (2.1e) 

y    yy 


x   =  w,  (2. If) 

xy     'xy 


where  e„  and  e„  are  the  extensional  strains,  v    is  the 
x      y  '  xy 

shear  strain,  Xv  and  X   are  the  curvature  changes,  Xv,7  is  the 
twist,  W  is  the  radial  displacement  of  the  midsurface  of 
the  nearly  circular  cylinder,  and  W  is  the  radial  deviation 
of  the  unstrained  midsurface  from  the  midsurface  of  a 
perfectly  circular  cylindrical  shell.   Hence 

W*  =  W  +  W 
Here,  and  in  the  sequel 

(  ,      =  ILL 


(  )     "  ^U- 

[     ' 'yy       2 

3y 


Rings 


The  rings  are  assumed  to  be  circular  in  shape  and 


r  • 

radially  thin  compared  to  their  radius  a  .   The  strain- 
displacement  relations  taken  for  the  rings  are  the  ones 
used  by  G.  D.  Galletly  [26]  in  his  investigation  of  the 
free  linear  vibration  of  ring-stiffened  shells.   The 
strains  of  the  centroid  of  the  ring  are  given  by 
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xl     r 
a 


r 
*2 


u 


♦    = 


(3  + 


(a1) 


r^2 


(-1 


(wr  +  w*ee) 


a 


(2.2) 


where  x-i  is  the  curvature  change  in  the  plane  normal  to  the 
ring,  x?  is  the  curvature  change  in  the  plane  of  the  ring, 

<}>  is  the  angle  of  twist  per  unit  circumferential  length, 

r 
and  £fl  is  the  circumferential  strain. 

Assuming  8  to  be  small,  e<<a,  and  that  W  and  W,   are 

zero  at  the  point  of  attachment,  the  displacements  and 

rotation  of  the  ring  can  be  geometrically  related  to  the 

shell  displacements  [26].   This  leads  to 


wr  =  W 


8  =  w, 


x 


y  (2.3) 


U   =  U  +  e 


Vr  =  V  +  e  W,. 


and  the  radius  of  the  ring  is  given  by 

ar  =  a  +  e 

Thus,  the  ring  strains  can  be  expressed  in  terms  of  the 

shell  displacements.   Substituting  (2.3)  into  (2.2)  and 

noting  that  (  ) ,   =  a(  ) ,   gives 

6        y 
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Xi  =  "4  [W,   +  S^  (U     +  e  W,    )] 
a  a 


X"  = 


<f>   = 


^2  (W  +  a2  W,yy) 


(a*) 


U,   +  e  W, 
—  (-W,    +  — X SSL) 


r  _   a 


w, 


eS  =  ^  (v'y  +  e  w'yy  "  ? 


f  (2.4) 


Note  that  (2.4)  does  not  contain  W,  the  initial  imper- 
fection of  the  shell.   This  is  due  to  the  fact  that  the 
intial  imperfection  W  and  its  slope  W,   have  been  taken 
equal  to  zero  where  the  rings  are  attached. 

IV.   CONSTITUTIVE  RELATIONS 

The  resultant  forces  and  moments  per  unit  length  of 

shell  midsurface  are  obtained  by  integrating  the  stresses 

over  the  thickness.   The  normal  forces  are  N   and  N  ,  the 

x      y 

shear  forces  are  N   and  N   ,  the  bending  moments  are  M 

xy      yx  *  x 

and  M  ,  and  the  twisting  moments  are  M„„  and  M„  .   The 
y  3  xy      yx 

positive  directions  of  these  forces  and  moments  and  the 
applied  pressure  P  are  as  shown  in  Figure  2.3. 

The  stress  resultants  are  related  to  the  strains  by 
the  relationships 


N  =  B(e   +  ve  ) 
x      x     y' 

Ny  =  B(ey  +  vex) 


N    =  N 
xy     yx 


_  1-v 


By 
2     XY 


(2.5a) 
(2.5b) 
(2.5c) 
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NXY  NYX 


FIGURE   2.3 


POSITIVE     DIRECTIONS    OF      FORCES    AND    MOMENTS 
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M   =  -D(X   +  VX  )  (2.5d) 

x       x     y 

My  =  "D(Xy  +  VXx)  (2.5e) 

M    =  -M    =  (1  -  v)  DX  (2.5f) 

xy     yx  Axy 

where  v  is  Poisson's  ratio,   and  B  and  D  are  the  extensional 
and  flexural  rigidities,  respectively,  defined  by 


B  =  -Jh 


1-v2 


D  =  Ehf 


12  (1-v2) 


M2.6) 


Young's  modulus  is  denoted  by  E. 

V.   EQUATIONS  OF  MOTION 

The  governing  equations  of  motion  are  derived  in  this 
section  by  means  of  Hamilton's  Principle.   For  the  case 
where  nonconservative  forces  are  present,  Hamilton's 
principle  states  that  [27] 

6  /  z    (3^+^K   =  0  (2.7) 

tl  dt 

where  3f°  is  the  total  work  function,  V  is  the  total  kinetic 
energy,  and  t  is  the  time  coordinate,  t,  and  t_  being  arbi- 
trary times. 

The  total  work  function  is  the  sum  of  the  work  function 
of  the  external  pressure  and  damping  force  9^?,  the  work 
function  of  the  shell  3^°  ,  and  the  work  function  of  the 


rings  3P  . 
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Thus 


0P=    W  +    W  +   3E  (2.8) 


The  variation  of  3/      is  given  by 

j-i 


L   2-rra 
6^°_  =  -  /   /    (C  W, .  -  P)  6W  dy  dx  (2.9) 

h  o   o        t 


where  C  is  the  viscous  damping  coefficient.   The  terms CW, 
and  P  are  held  constant  during  the  variation. 
The  variation  of  9f°    is  given  by 

L      2-rra 
s  oo         Lxx  xy'xy  yy 


(2.10) 


-M^,6y   +2M   6y    -  M  6y    dydx 
x  Ax     xy  Axy    Y     Y 


The  in-plane  forces  and  moments  are  also  held  constant 
during  the  variation. 

The  contribution  Sir     of  the  ring  forces  to  the  total 
work  function  is  taken  as  the  negative  of  the  strain  energy 
in  the  rings.   Hence,  6df     is  given  by  [26] 

Nr    27Tar 


1  =  1   ° 


2  r-  2 

+  G  C^  (<j>)   +  E  A   (ej)  ]  dy  6(x-x.)  dx 
r  i  -L  r    y  ^ 


(2.11) 
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where 

y       characteristic  length 

Iv,  I   moments  of  inertia  of  the  ring  cross  section 
x    z 

with  respect  to  axes  passing  through  the 

centroid  of  the  ring 
A       cross  sectional  area  of  the  ring 
CT      St,  Venant ' s  torsion  constant  for  the  ring  [28] 
6(  )     Dirac  delta  function,  not  to  be  confused  with 

the  symbol  used  to  denote  a  variation 

N       number  of  rings 
r  ^ 

E       Young's  modulus  of  the  ring  material 
Gr      shear  modulus  of  the  ring  material 
Pr      mass  density  of  the  ring  material 
j        subscript  which  associates  a  variable  with 
the  j    ring;  for  example,  x.  is  the 

x-coordinate  of  the  j    ring 

—   =r 

dy  =  fL_  dy 
a 

The  Dirac  delta  function  has  been  introduced  here  to  allow 

the  integration  of  69f°  to  be  carried  out  over  the  length 

of  the  shell.   This  technique  has  been  used  by  Block  [23] 

and  is  necessary  in  order  to  combine  3^?,  2HQ,    and  Z/fL    in 

a  double  integration  form.   It  requires  the  introduction 

of  a  characteristic  length  y  that  corresponds  to  the  width 

of  a  ring. 

The  total  kinetic  energy  of  the  structure  is  given  by 

=  Sri   +  &.  (2.12) 


V  =  7  s   +  ^ 
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where  %f      is  the  kinetic  energy  of  the  shell  and  %f    is  the 
kinetic  energy  of  the  ring.   Only  the  radial  velocities  will 
be  included  in  the  kinetic  energy  of  the  shell  and  of  the 
ring;  the  axial,  circumferential,  and  rotatory  inertia  of 
both  the  rings  and  the  shell  are  assumed  to  be  negligible. 
Thus,  the  kinetic  energy  of  the  shell  is  given  by 


nh   L   27ra     o 

&   * =^p  /   /    (W,.  r  dydx  (2.13) 

v  s    ^  o  o      t 


where  p  is  the  mass  density  of  the  shell  material.   The 
kinetic  energy  of  the  rings  is  given  by 

.         PrAr     Nr       L      27rar  2 

^C   =   ~9Tr      Z      ;      f  <W'J       d^    Mx-x.)    dx  (2.14) 

*  ^M      j  =  i   o      o  r  J 

Substituting  (2.8)  and  (2.12)  into  (2.7)  yields 

t2 
6/   (3#  +  WL   +     9P   +  g*   +  Wi)    dt  =  0        (2.15) 

•f-         J-r         S         T  o         L 

Ul 

Substituting  (2.1)  into  (2.10);  (2.4)  into  (2.11);  and 
(2.9),  (2.10),  (2.11),  (2.13),  and  (2.14)  into  (2.15)  and 
performing  the  variation  yield  the  following  equations  of 
dynamic  equilibrium 

N 
r  EI 

N     +  N      +   Z  [-r- -  (W,VT„r  +  aU, 

x,x   xy,y   .  1l  ay    '*yy     yyyy 

G  C  (2*16) 

+  ea  W,      )  .  -  -f~  (i  U,    -  W,    )  .]6(x-x  )  =  0 
xyyyy  j     ay    a    yy      xyy  -j       j 
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Nr  e   A 

Nv   __    +   N         v    -    I   -E-£    (V.         -  e  W, 
y,y  xy,x      ,    ^    y  yy  yyy 


-   -  W,    ) .    6(x-x.)    =   0 
ay:  3 


(2.17) 


DV4W   =    N„(W,  +    W,       )     +    2N    „     (W,  +    m,       ) 

x v    ' xx  'xx'  xy         'xy  xy 

Nrp   A 
+   Ny(W,H+W,yy)+-Ny  +  P-   CW,t    -phW,tfc   -   ^-^(w.tt). 

Nr     E    I  _  , 

^y^T    <w'xx  +  aU'xyy  +   2ea  w'xxyy  "  ea     "'xyyyy 

—o    o  ^r^x         9  4 

+  ezaz   W,       )  . j-  (W  +  2az  W-  „7  +  a  W,     )  . 

xxyyyy'j    ya4  'yy        YYYY   J 

E  A  2 

+  *    (ea2  V     +  a  V,„„  -  ea  W,     -  2ea  W/vv-W) . 

2      /yyy     yy       yyyy       W  '  j 

y  a 

GrCT 
+  -77T-  (U,xvv  -aW,     )  .]  6(x  -  x.) 

ya      XYY        xxyy  3  D  (2.18) 


VI.   BOUNDARY  CONDITIONS 

The  variational  procedure  also  yields  the  admissible 
boundary  conditions  [29] . 

Simply  Supported 

The  boundary  conditions  that  correspond  to  simply 
supported  ends  are 
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V 

=  0 

w 

=  0 

(2 

M 

X 

=  0 

N 

X 

=  0 

at  x 

=i 

0 

,L 

Clamped 

The 

boundary 

cond 

itions 

that  < 

correspond 

to 

clamped 

ends  are 

U 

=  0 

V 

=  0 

(2 

W 

=  0 

w, 

X 

=  0 

at  x 

= 

0 

,L 

VII.   NONDIMENSIONAL  EQUATIONS 


Nondimensional  Variables 


The  nondimensional  variables  and  parameters  that  will 
be  used  in  the  analysis  are  defined  by 

\ 


e-f 


t      F     % 
T  =  -  ( Z—    ) 

a  v  ,,   2  ' 


L 
a 


P(l-v") 


u 

u  =  a 


p  =  F*d-V  > 


V 
v  =  a 


Eh 


« ■ c  ^h 


=   Y. 
a 


W 

w  =  I 


-        W 

w  =  — 

a 


(2.21a) 
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t,     X 

i   2 

Eh     ne  • 

=  Ny 

1-v2 
Eh 

n?e  - 

i   2 
xy  Eh 

h 
a  =  a 

e 

e  =  — 
a 

1     E 

Ix  (1-v2) 
a  hy 

Pr 
2     p 

Ar 
hy 

C   =  ^ 
3    E 

Ar  (1-v2) 
hy 

Er 

C   = 

I   (1-v2) 
z 

4     E 

2u 

a  hy 

5     E 

CT  (1-v2) 

i„2  u,. 

(2.21b) 


Nondimensional  Equations  of  Motion 

Expressing  the  equations  of  motion  (2.16),  (2.17), 
and  (2.18)  in  terms  of  the  above  non-dimensional  quantities 
leads  to  the  three  equations 


Nr 

V?  +  nce,e  \l=1[c4    (w'cee  +  u'eeee  +  ^'ceeee^ 


(2.22a) 

J 


nQ0   +   nro    r         S    C,     (w,_    +   ew,nAn    -    v,flfl).    5  (£-£  . )    =    0 


»,e  '   iAce,c   '^ 


(2.22b) 
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12      4  -  - 

12    a      v      w   =   n^     (w,^+w,^)    +    2n^e     (w,^e+w,^0) 


N. 


c    . 


+    nfl     (w,fifi+w,Afi)    +    nft    +   p  -w   -[w+    Z    C7   w^    6 (£-£  J  ] 


a    "       1"'J:1V'2      : 


N, 


+    Z     [C,     (w,OQ+u,         +2ew,  +eu, 


j  =  l 


C^ee       '?' 


(2.22c) 


+   e      w,__         )  .    -    C,     (w   +    2    w,         +  w,  ) 

ueeee  d        i  'ee  eeee  j 


+    C-,     (ev,  +   v.       -    e      w,  -    2e   w,         -w) 

3  ftflfl  ft  flftflfl  oft 


+    Cc     (u,  -   w,  )     ]     6  (£-£  .  ) 


Here    (')    =    8 (    )/9t. 

The  in-plane    forces    N.N.    and   N„.,    can   be    expressed    as    func- 

x      y  x^ 

tions  of  the  displacements  U,  V,  and  W  by  combining  the 
constitutive  relations,  (2 . 5a) -  (2 . 5c)  ,  and  the  strain- 
displacement  relations  (2 . la) -  (2 . lc) .  In  nondimensional 
form 


ru  = 


n„  = 


n. 


2    -  2    _ 

U,£  +  h(W,r)  +     W,rW,r     +V [V,Q-W+^ (W,Q)    +  W , QW , Q ] 

(2.23a) 

2    —  2  — 

v/Q  -  w  +  5§(w,e)   +  w,6w,Q  +  v[u,^+1^(w/  )  +W,  W,  ] 


(2.23b) 


1-v  .  - 

~Y~    (v,_  +  u,Q  +  w,  w,Q  +  w,  w,Q  +  w,Qw,  ) 


(2.23c) 
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Nondimensional  Boundary  Conditions 

Simply  Supported.   The  boundary  conditions  for  the 
simply  supported  shell  are  expressed  by  (2.19).   The 
corresponding  conditions  on  the  nondimensional  displace- 
ments v  and  w  are 

v  =  0  (2.24) 

w  =  0  at  £  =   0,1  (2.25) 


The  bending  moment  M   is  a  function  of  the  curvatures  and 
was  expressed  by  (2.5e)  as 


Mx  =  "D(xx  +  VV  (2.5e) 

Expressing  the  curvatures  in  terms  of  the  displacements  by 
using  (2. Id)  and  (2.1e)  gives 


M   =  -D(W.  „  +  vW,   ) 

x       xx     yy 

At  the  ends  of  the  shell  W  =  0  which  implies  that  W,    =0; 

yy 

hence,  the  condition  that  M  =  0  at  x  =  0  and  x  =  L  implies 

x 

that  W,    =  0  at  x  =  0  and  x  =  L.   In  nondimensional  form 
xx 

w,    =0  at  ^  =  0,i  (2.26) 

Since  N  =  0  at  x  =  0  and  x  =  L,  it  follows  that  n   =  0  at 
x  £ 

£   =    0    and  £  =  i.   Hence,  from  (2.23a) 

2    —  2   — 

n£  =  U'£  +  %(™'e)       +   w'rw'r  +  v[v,  -w+?g(w,  )   +  w,Qw,Q]  =  0 

atC=  0,i  (2.27) 
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Clamped.   The  clamped  boundary  conditions  (2.20)  can 
also  be  expressed  nondimensionally  as 


u 

=    0 

V 

=    0 

w 

=    0 

w,, 

-    =    0 

at  £  =  0,i 


(2.28) 
(2.29) 
(2.30) 
(2.31) 


Response  Symmetric  to  Midspan.   When  the  load ,  the 
initial  imperfections,  and  the  response  are  symmetric  with 
respect  to  the  midspan  of  the  shell,  only  the  portion  of  the 
shell  between  £  =  0  and  £  =  i/2  needs  to  be  considered.   The 
nondimensional  conditions  that  express  symmetry  at  the 
midspan  are 


u 

v, 


w 


'K 


w, 


KU 


at  K    =    x/2 


(2.32) 
(2.33) 
(2.34) 
(2.35) 


The  three  equations  of  motion,  (2.22),  together  with 
the  expressions  for  nr,    na,    and  n_Q,  (2.23),  and  the  boun- 
dary  conditions  constitute  a  formulation  of  the  problem 
under  investigation.   The  next  two  chapters  are  devoted  to 
their  solution. 
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CHAPTER  III 
DEVELOPMENT  OF  THE  NUMERICAL  MODEL 
I .   INTRODUCTION 

In  this  chapter  the  nondimensional ,  nonlinear  govern- 
ing partial  differential  equations  derived  in  Chapter  II 
are  replaced  by  a  set  of  difference-differential  equations 
in  which  time  is  the  only  remaining  independent  variable. 
The  difference-differential  equations  are  henceforth 
referred  to  as  the  numerical  model. 

The  model  is  developed  by  expanding  the  dependent 
variables  in  Fourier  series  in  the  circumferential  coordi- 
nate.  The  Fourier  coefficients  in  the  series  are  functions 
of  t  and  the  axial  coordinate  £.   Products  of  Fourier  series 
arise  wherever  there  are  nonlinear  terms  in  the  governing 
equations.   The  expansion  of  these  products  to  single  series 
generates  products  of  the  Fourier  coefficients  that  couple 
the  modes.   The  derivatives  of  the  Fourier  coefficients 
with  respect  to  the  axial  coordinate    are  then  approximated 
by  modified  finite  differences,  leading  to  the  set  of 
difference-differential  equations . 

II.   FOURIER  SERIES  EXPANSIONS 

During  the  early  phases  of  the  present  study,  con- 
sideration was  given  to  a  numerical  model  that  approximated 
both  the  circumferential  derivatives  and  the  axial 
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derivatives  by  means  of  finite  differences.   To  achieve  good 
resolution  of  the  strains  and  curvatures  of  the  cylinder, 
the  grid  points  must  be  close  enough  together  so  that  there 
are  several  grid  points  covering  the  shortest  wave  length 
of  the  response.   Based  on  the  works  of  previous  investi- 
gators [9,  10,  20],  the  parametrically  excited  modes  are 
expected  to  have  many  waves  around  the  circumference. 
Consequently,  a  very  large  number  of  points  is  needed  for 
adequate  resolution  of  the  parametrically  excited  modes. 
For  example,  a  flexural  mode  with  sixteen  waves  around 
the  circumference  is  not  unreasonable  to  expect.   Taking 
the  number  of  axial  stations  to  be  twenty  and  assuming  that 
eight  node  points  per  wave  length  are  needed  around  the 
circumference,  the  total  number  of  node  points  is  approxi- 
mately 2500.   Satisfying  axial,  circumferential,  and  radial 
equilibrium  at  each  node  yields  over  7500  simultaneous 
equations.   The  number  of  equations  can  be  reduced  by  one 
half  if  the  response  is  assumed  to  be  symmetric  with  respect 
to  some  axial  plane.   Even  so,  the  number  of  equations  is 
still  very  large. 

However,  if  the  displacements  are  expanded  in  Fourier 
series  in  the  circumferential  coordinate,  instead  of  using 
finite  differences,  the  circumferential  variation  of  the 
normal  modes  of  the  cylinder  can  be  represented  exactly. 
Furthermore,  those  modes  which  do  not  respond  can  be  ex- 
cluded. Thus,  the  number  of  equations  can  be  kept  small 
while  retaining  the  significant  features  of  the  coupling 


58 


between  modes.   Expansions  of  the  type  employed  here  have 
been  used  previously  by  R.  E.  Ball  [30]  in  a  nonlinear 
analysis  of  statically  loaded  shells  of  revolution. 

Form  of  the  Series  Expansions 

The  assumption  is  made  that  the  initial  imperfection 
of  the  cylinder  can  be  expressed  by  the  Fourier  cosine 
series 

00 

w(£,6)  =   E  wn  (K)    cos  n9  (3.1) 

n=0 

This  form  somewhat  limits  the  applicability  of  the  analysis 
since  the  initial  imperfection  of  most  shells  cannot  be 
expressed  by  the  cosine  series  alone.   However,  the  form 
is  still  sufficiently  general  to  permit  the  retention  of 
coupling  effects  that  hitherto  have  not  been  retained. 
There  are  no  conceptual  problems  involved  in  using  the 
complete  Fourier  series.   Only  the  cosine  series  is  used 
here  because  of  limitations  imposed  by  the  practical  con- 
siderations of  computer  storage  space  and  computing  time. 
In  order  to  correspond  with  the  form  of  w  given  by 
(3.1),  the  remaining  dependent  variables  are  expressed  in 
trigonometric  series  as  follows: 

00 

u(£,9,t)  =   Z   un  (£,t)  cos  n6  (3.2a) 

n=0 

oo 

v(£,8,t)  =  Z      vn  (£,t)  sin  n6  (3.2b) 

n=l 

oo 

w(£,G,t)  =   £   wn  (£/T)  cos  n6  (3.2c) 

n=0 
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n  (£,e,x)  =   Z   n^  (£,t)  cos  n6  (3.3a) 

5  n=0   ^ 


ne(£,0,T)  =   Z   nQ  (£,t)  cos  n0  (3.3b) 

n=0 

00 

n   (€,6,t)  =   Z   n"  (£,t)  sin  n6  (3.3c) 

t.0  n=l   ^e 


When  these  expressions  for  the  displacements  and  in- 
plane  forces  are  substituted  into  (2.22)  and  (2.23),  the 
linear  terms  can  be  grouped  together  in  one  series.   However, 
the  nonlinear  terms  give  rise  to  products  of  series  and 
cannot  be  grouped  with  the  linear  terms  until  they  are 
expanded  into  single  series.   The  products  of  series  are 
of  three  types;  products  of  (1)  two  cosine  series,  (2)  a 
sine  series  and  a  cosine  series,  and  (3)  two  sine  series. 
The  details  of  the  expansion  of  these  types  of  products 
are  described  in  Appendix  B  of  Reference  [30].   Employing 
(3.1),  (3.2),  and  (3.3),  the  nonlinear  terms  contained  in 
(2.22c)  and  (2.23)  can  be  expressed  by 


2  "i  n  °°   n 

^(w,^)   =  \  (  £  w,r  cos  i9)  (  Z  w,^cos  j0)=  Z  C,rv  cos  n9 

(3.4a) 


i=0   ^         j=0   ^       n=0  xx 


oo 


2  i  -i  n 

^(w,  )   =  H  (  Z  -iw  sin  i6)(  Z  -  jwJsin  j0)  =  Z  CTm  cos  n9 

6        i=l  j=l  n=0 

(3.4b) 


w,,w,fl  =  (  Z  w)  cos  i0) (  Z  -jw^sin  j0)=  Z  CXT  sin  n0 
^        i=0   ^         j=l  n=l 

(3.4c) 
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00  00  00 

w,    w,    =(    Z   w,J    cos    i0) (    Z   w,3    cos    j6)=    Z   C  cos   n6       (3.4d) 

**      %      i=0      *  j=0      *  n=0    IXX 

00  00  00 

w,0w,0=(    Z    -iw1    sin    i0)  (    E    -jw-3    sin    j0)=    Z    CTrTur   cos   n0 
6      6      i=l  j=l  n=0    iri 

(3.4e) 

00  00  00 

w,rw,0=(    Z   w, r   cos    i0) (    Z    -jw      sin    j0)=    Z    Cj„T    sin   n0 
i=0  j=l  n=l 


(3.4f) 


(3.4g) 


(3.4h) 


00  OO                                                                               00 

■j  -i 

w,    w,    =(    Z    -iw  sin    i0) (    Z   w,  cos    j0)=    Z    CT_V   sin   n0 

s      1=1  j=0      s                       n=l 


OO  00 

nr (w,rr+w,rr)=  (Z   n1   cos    i0) [    Z    (w,J    +w,Jr)cos    j0] 
^        ^        ^      i=0    *  j=0        "        ^ 


=    Z    rir    cos    n( 

n=0    * 


00  OO 

2n       (w,^e+w/^e)  =  (    Z    2n^e    sin   i0)[    Z    -j  (w,£+w,jj)    sin    j0] 

(3.4i) 

00 

=     Z     Tiro     COS     n0 

n=0    C6 


00  o 

--1        fw^4.^ 


nfi(w'flo+w'flJ  =  (    z    nfl    cos    ie)  t    Z    -j     (wJ+wJ)    COS    j0] 
o         00         00         i=Q  j=Q 


00 

v         n  o 

=      Z    n_    cos   n0 
n=0    6 


(3.4j) 


The  products  containing  the  initial  imperfection,  or  deri- 
vates  of  the  initial  imperfection,  are  not  nonlinear  since 
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w(£,6)  is  assumed  to  be  known.   However,  the  products  con- 
taining w  expand  in  the  same  way  as  the  nonlinear  products; 
hence,  they  are  included  in  the  above  list. 

Reduction  to  Two  Independent  Variables 

Equations  of  Motion.   Expressions  for  the  Fourier 
coefficients  of  the  in-plane  forces  n^  ,    nj},  and  n.fl/  can  be 
obtained  by  substituting  (3.2),  (3.3)  and  (3 . 4a) - (3 . 4g)  into 
(2.33),  and  equating  coefficients  of  like  values  of  n. 
Performing  the  operations  yields 


n*  =  u  "  +  v(nvn  -  wn)  +  c£x  +  C£xx 

v (C^T  +  C5TT)         n  =  0,  1,  2,  .  .  . 


n*  =  nv11  -  wn  +  vu  "  +  c£m  +  Cn 
6  £     TT     ITT 


+    V(CXX   +   C?XX)  n    =    0,    1,    2,     . 


(3.5a) 


(3.5b) 


n  1-v 

n 


(v   n    -    nu11   +   Cn      +   Cn         +   C11       ) 
£6  2     (V'^         nu      +    LXT    +    LIXT         LITX; 

n   =    1 ,    2 ,    ... 


(3.5c) 


The  equations  of  motion  in  the  axial  and  circum- 
ferential directions  can  be  expressed  in  terms  of  the 
Fourier  coefficients  by  substituting  (3.2)  and  (3.3)  into 
(2.22a)  and  (2.22b),  and  by  making  use  of  (3.5).   The 
result  is 


62 


u  *      +±±Z     nv"   -i=l     n2un   +    Xln   +      S    [CAn4un   + 
55        2  5        2  j=1      4 

C4n2    (en2-l)    wr°   -   C5n2    (w,£   -   un) ]  §    6(g-£j)    =   0 


n   =    0 ,    1,2,    ...  (3.6) 


2    n         1+v   _      n         1-v        n       ,     % 0n 
5  2      V'SS 


-  nu/c    +   — —  v,rr    +    A2 


Nr 

+      E      CJn(en   -1)    wn   +   n2   vn]       6  (£-£    )    =    0 
j=l      3  J  J 

n=l,2,     ...  (3.7) 

where 


(3.8) 


and 


+  — —  n(Cn      +   Cn         +   Cn       i    -    vw11 
2      n^XT  IXT         UITX;  VW,C 


A2n   =    -n[Cn      +   Cn         +    v(Cn      +   C^vv)] 
TT  ITT  XX  IXX' J 


+    1-v    tnn  n  n  n 

+   "T"    (CXT,C         CIXT,5         CITX,C)  nW 


(3.9) 


All  of  the  nonlinear  terms  are  contained  in  XI   and  X2    . 
Examination  of  these  two  quantities  reveals  that  the  axial 
and  circumferential  equations  are  nonlinear  in  w  only; 
there  are  no  nonlinear  terms  which  contain  u   or  vn. 

By  substituting  (3.2c),  (3.3b),  and  (3 . 4h) -  (3 .4j )  into 
(2.22c),  the  equation  of  radial  equilibrium  is  obtained  in 
the  form 
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1   2  ,   n      n    2  4   n  r..n 

T2  a   (W'^£  "  2n   W'?C  +  n   w  )  =  6Qn  p  -  [w   + 

Nr 

Z   C_  w?  6  (£-£.)]  "  |  wn  +  n£ 
j  =  1   2   :       D      a       « 

Nr 

+  n?  +  n?fl  +  n"  +  z  r"  6(C-5.)  (3.10) 

n  =  0  ,  1,  2 ,    ... 
where 

Rj  =  -[C4  n2  +  C1  (n2-l)2  +  C3  (en2-l)2]w!? 


-[C,n2(l-en2)  +  Crn2](u?r)   -  C  n  (en2-l)vn      (3.11) 
4  J       s,  j     3  j 

-[C.en2(2-en2)  -  C  ch2] (w*   ) . 

and 


1  if  n  =  0 
un    0  if  n  ?   0 


The  number  of  independent  variables  has  now  been  re- 
duced from  three  (£,  0,  and  t)  to  two  (£  and  t).   However, 
now  there  are  three  equations  to  be  solved  for  each  value 
of  n  considered. 

Boundary  Conditions 

The  boundary  conditions  must  also  be  expressed  in 
terms  of  the  Fourier  coefficients  of  the  displacements. 
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Simply  Supported.   Substitution  of  (3 . 2b) -  (3 . 3a)  into 
(2 . 24 ) -  (2 . 27)  yields  the  boundary  conditions  for  a  simply 
supported  shell 

vn  =0  (3.12a) 

wn  =0  (3.12b) 

w?  =  0  (3.12c) 

n"  =0  at  E,    =    0,i 

With  (3.5a),  (3.12a),  and  (3.12b)  the  last  of  the  above 
conditions  reduces  to 

[u,£  +  c£x  +  c£xx  +  v(c^T  +  c^TT)]  =0     at  K  =  0,1 

From  equations  (3.12b),  (3.4b),  and  (3.4e),  it  follows  that 
C^„   and  Cj        are  zero  at  the  ends  of  the  shell.   Hence, 
the  final  boundary  condition  becomes 

[u'?  +  CXX  +  CIXX  ]  =  °        at  C  =  0,1        (3.12d) 

Clamped.   The  boundary  conditions  for  the  clamped 
shell  are  obtained  by  substitution  of  (3.2)  into  (2.28)- 
(2.31) .   The  result  is 


un   =  0  (3.13a) 

vn   =  0  (3.13b) 

wn   =  0  (3.13c) 

wn.  =0  at  £  =  0,1         (3.13d) 
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un 

=    0 

n 
V    '5 

=    0 

wn,^ 

=    0 

wn,,, 

=    0 

Response  Symmetric  to  the  Midspan.   Substitution  of 
(3.2)  into  (2.32)- (2.35)  gives 

(3.14a) 

(3.14b) 

(3.14c) 

,£     -v  at  £  =  0,  i/2      (3.14d) 

III.   MODIFIED  FINITE  DIFFERENCES 

In  this  section  a  brief  discription  of  the  differenc- 
ing procedure  used  to  approximate  the  axial  derivatives  of 
(3.6)- (3.14)  is  presented.   These  finite  differences  have 
been  called  "modified"  differences  [31,  32]  or  "half 
station"  differences  [33] . 

As  an  alternative  to  the  finite  differences,  the  de- 
pendent variables  could  have  been  expanded  in  Fourier  series 
in  the  axial  coordinate.   T.  Wah  and  W.  C.  L.  Hu  [34]  made 
use  of  expansions  of  this  type  in  their  linear  vibration 
analysis  of  simply  supported,  ring-stiffened  cylinders. 
The  equations  of  motion  of  the  shell  and  the  equations  of 
motion  of  the  rings  were  satisfied  simultaneously  by  enforc- 
ing the  continuity  of  displacements,  forces,  and  moments 
between  the  shell  and  the  rings.   However,  the  present 
problem  is  more  complex  than  their  problem  due  to  nonlin- 
earities  in  the  equations  of  motion  of  the  shell  and  due  to 
the  clamped  boundary  conditions.   The  use  of  the  finite 
differences  to  approximate  the  axial  derivatives  appears 
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to  facilitate  the  handling  of  both  the  boundary  conditions 
and  the  effects  of  the  discrete  ring  stiffeners.   The  use 
of  series  expansions  in  the  axial  coordinate  was  not  in- 
vestigated in  detail;  consequently,  the  advantages  and 
disadvantages  are  not  clear  at  the  present  time. 

Modified  differences  are  applicable  to  systems  of 
equations  in  which  only  even  or  odd  ordered  derivatives 
of  a  given  dependent  variable  appear  in  any  one  equation. 
As  an  example  of  equations  that  fall  into  this  category, 

consider  (3.6)- (3.10)  with  nn  replaced  by  (3.5b)  and 

6 

without  the  nonlinear  terms  and  the  terms  associated  with 
the  rings ,  that  is 

U^C  +  ^r"  nV^r  "  HF  n2  ^  "  V  W^  =  °  (3.15a) 

n  =  0 ,  1 ,  ... 

2  n    1+v     n     1-v   n        n 
-n  v   -  — —  n  u,>.  +  — —  v,    +  n  w   =0         (3.15b) 

n  =  1 ,  2 ,  ... 

17  a   (W'EEE£  ~  2n   W'rf  +  n   w  )  (3.15c) 

r        -n    c  »n     n    n     n 

=  o    p-w   -  —  w   +  nv   -w   +  vu,r 
On  a  '£ 

n  =  0 ,  1 ,  ... 


Considering  the  spatial  derivatives,  note  that  only 
even  derivatives  of  u   appear  in  equation  (3.15a);  only  odd 
derivatives  of  un  appear  in  (3.15b)  and  (3.15c).   Similarly, 
equation  (3.15a)  contains  only  odd  derivatives  of  v   and  w  ; 
and  (3.15b)  and  (3.15c)  contain  only  even  derivatives  of 
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n      n 
v   and  w  .   For  systems  of  equations  having  this  property, 

a  staggered  finite  difference  mesh  can  be  used  advan- 
tageously.  A  staggered  mesh,  as  illustrated  in  Figure  3.1, 
is  one  where  half  stations  are  interspaced  between  the 
whole  stations,  and  the  distance  between  consecutive  whole 
or  half  stations  is  denoted  by  d.   In  the  modified 
difference  method,  some  of  the  dependent  variables  are 
defined  only  at  whole  stations  and  the  remainder  are  defined 
only  at  half  stations.   For  reasons  which  are  not  yet 
obvious,  u  will  be  defined  at  half  stations;  while  v   and 
w   will  be  defined  at  whole  stations  as  shown  in  Figure  3.1. 


1 
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Figure  3.1 
STAGGERED  FINITE  DIFFERENCE  MESH 
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K.  P.  Chuang  and  A.  S.  Veletsos  [31]  were  apparently 
the  first  to  employ  the  modified  differencing  method  in 
analyzing  circular  cylindrical  shells.   They  proposed  the 
fundamental  modified  central  difference  operator 
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[(  ),  1   =  i  [(  ),    -  (  )'   J  (3.16a) 

for  derivatives  with  respect  to  the  independent  variable  x. 
Consistent  applications  of  the  operator  yield   a  set  of 
difference  approximations  for  the  higher  order  derivatives 
of  the  form 

(f,   ),  -  -hr   (f    -  2f.+  f    )  (3.16b) 

'xx  k   d2    k+1     k    k-l 

(f,xxx),  =  -\    (f..3  "  3f     +  3f     -  f   3)     (3.16c) 
xxx  k   d3   k+-     k+i^     k-Jg    k-J 

<f'xxxx>k  =  Jr    <fk+2    -    "k+1   +   6fk   -    4fk-l   +   fk-2> 

(3.16d) 

Returning  now  to  (3.15),  the  difference  equations  are 
obtained  by  expressing  (3.15a)  at  half  stations  and  (3.15b) 
and  (3.15c)  at  whole  stations.   The  resulting  difference 
equations  are 

|i  K+f  -  HU  +  U!U>  +  ¥  4  <vx+rvx)      b-w 
-  t1"2  uk+% "  vi  (wk+i  -  wk'  ■  °    n  =  °'x---- 

-n2        v"  -    liS  n   i    (u"  .    -u?  .  )+l=i    1     <v"      -2v?+v"    ,) 

k  2  d        k+^        k-V       2      d7  k_1 

+   nw£   =    0  n   =    1,    2,    .  .  .     (3.17b) 


12rl.n  n  n  _  n      .       „   2    1  ,   n        ^   n 

12   a    [^4    (wk+2    "    4wk+l   +    6wk   "    4wk-l+  Wk-2>-    2n      ^2(Wk+r2"k 

x    n      \     .      4   n,         ~  ..n        c    *n  n  n  1,   n  n     \ 

+  wk_l)    +  n  wk]    =    6^  pk   -  wk   -  -  wk   +  nvk   -  wk   -f   vjfu^-u,^) 

n   =    0,    1,    ...  (3.17c) 
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Observe  that  in  the  development  of  these  three  equations, 
it  was  not  necessary  to  define  un  at  whole  stations.   Like- 
wise, w  and  v  did  not  require  definition  at  half  stations. 
These  circumstances  prevail  for  two  reasons:  (1)  the  deriv- 
atives of  any  one  variable  are  all  odd  or  even  ordered 
within  a  given  equation,  and  (2)  the  equations  are  linear. 

Superiority  of  Modified  Differences  over  Conventional 
Differences 

Chuang  and  Veletsos  [31]  showed  that  when  the  con- 
ventional central  differences  are  used  for  all  the  deriv- 
atives, the  difference  equations  obtained  by  discretization 
of  the  governing  differential  equations  do  not  agree  with 
those  obtained  by  minimizing  the  finite  difference  form 
of  the  potential  energy.   On  the  other  hand,  when  modified 
differences  are  used,  the  difference  equations  obtained  by 
both  procedures  are  identical.   Consequently,  the  conven- 
tional and  modified  differencing  of  the  differential  equa- 
tions gives  rise  to  two  different  numerical  models,  a 
conventional  difference  model  that  approximates  only  the 
differential  equations,  and  a  modified  difference  model 
that  approximately  satisfies  both  the  differential  equa- 
tions and  the  condition  of  minimum  potential  energy.   In 
Appendix  A,  this  inconsistency  between  the  two  schemes  is 
shown  to  prevail  even  after  the  displacements  have  been 
expanded  in  Fourier  series  in  the  circumferential  coordinate 
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Chuang  and  Veletsos  [31]  compared  the  accuracy  of 
the  modified  and  conventional  differences  when  used  in 
the  static  analysis  of  cylindrical  shell  roofs.   They 
found  that  the  modified  differences  were  far  more  accurate 
than  the  conventional  differences  for  the  same  mesh 
spacing.   Ball  [32]  used  both  the  conventional  and  modified 
differences  to  compute  the  natural  frequencies  of  circular 
rings  and  compared  the  numerically  computed  frequencies 
with  the  exact  frequencies.   The  frequencies  obtained 
using  modified  differences  were  considerably  more  accurate 
than  those  obtained  with  conventional  differences  for  the 
same  mesh  spacing.   N.  J.  Cyrus  and  R.  E.  Fulton  [33] 
developed  analytical  representations  of  the  errors  asso- 
ciated with  the  modified  and  conventional  differences 
when  applied  to  beam  and  string  problems.   For  a  fixed 
number  of  stations,  the  errors  in  calculated  deflections 
and  curvatures  resulting  from  the  use  of  modified  differ- 
ences were,  in  general,  smaller  than  those  obtained  by 
using  conventional  differences. 

Adaptation  to  the  Nonlinear  Problem 

The  discussion  of  the  modified  differences  has  avoided 
applications  to  nonlinear  problems  until  now. 

Finite  Difference  Grid  and  Notation.   In  adapting 
the  modified  differences  to  the  nonlinear  equations,  (3.6)- 
(3.14),  the  pattern  used  for  the  linear  equations  is 
followed.   Difference  equations  are  written  which  make  use 
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of  a  staggered  node  system,  and  the  displacements  are 
defined  at  either  half  or  whole  stations.   A  convention 
is  adopted  whereby  both  half  and  whole  station  quantities 
are  given  integer  subscripts.   To  distinguish  whole  station 
quantities  from  half  station  quantities,  an  upper  case 
letter  will  be  used  for  the  whole  station  subscript,  and 
a  lower  case  letter,  for  the  half  station  quantities. 
As  shown  in  Figure  3.2,  the  k    half  station  node  is  located 
between  whole  station  nodes  K  and  K+l.   The  distance  between 
adjacent  whole  station  nodes,  as  well  as  the  distance  be- 
tween adjacent  half  station  nodes,  is  equal  to  d.   Let  the 
radial  and  circumferential  displacements  w   and  v   be 
defined  at  the  whole  stations;  and,  the  axial  displacement 
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Figure  3.2 
SPATIAL  DISCRETIZATION  USING  MODIFIED  DIFFERENCES 
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u  ,  at  half  stations.   At  whole  stations  the  difference 
operator  of  the  modified  differences  becomes 

[(  »'C]K=(i  [<  >k"  (  >k-d 
and  at  half  stations  it  becomes 

I(  »'?]k=d  [<  >K+1  "  (  'k1 

Here,  and  in  sequel,  k  =  K  within  a  difference  expression 
or  equation. 

Difference  Expressions  for  Derivatives.   Since  the 
three  displacements  u  ,  v  ,  and  wn  are  not  defined  at  every 
node  point,  the  derivatives  can  be  evaluated  at  either 
half  stations  or  whole  stations,  but  not  at  both.   The 
derivatives  that  can  be  evaluated  at  the  half  stations  are 

<VVk  =  I  (VK+1  "  VK»  (3-18a) 

<W'5>k    "   I    <WK+1    "    WK>  (3-18b> 

<V5'k  =  £  K+i  -  K  +  Vi>  <3-18c) 

'"'SCc'k  =  P  (WK+2  "  3WK+1  +  K    -    WK-1>        <3-18d» 

The  derivatives  that  can  be  evaluated  at  the  whole  stations 
are 

'"'s'k  =  I   (uk+l  "  uv»  (3-18e) 
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(wlWK    =    72     (WK+1    "    2wK    +   WK-1}  (3.18f) 

d 

(VVk   ■    3    (VK+1    "    2VK    +   VK-l'  <3-189) 


(w1}     )   =  —  (wn  n    -    4wn    +  6wn  -  4wn  ,  +  wn   ) 
1  'U^  K    d4  V  K+2    *WK+1      K    *WK-1    WR-2J 

(3.18h) 

Evaluation  of  Nonlinear  Terms.   The  nonlinearities 

in  the  axial  and  circumferential  equations  (3.6)  and  (3.7) 

are  embodied  in  the  functions  Al   and  A  2  ,  given  by  (3.8) 

and  (3.9);  in  the  radial  equation  (3.10)  the  terms  n^,  n*?, 

0    6 

n  r\ 

rirQ,  and  ne  are  nonlinear.   The  computation  of  these  non- 
linear quantities  in  terms  of  the  discrete  displacements 
presents  difficulties  not  encountered  with  the  linear 
equations . 

To  demonstrate  the  difficulty,  assume  that  A2   must 
be  computed  at  the  whole  station  K.   From  (3.9),  note  that 
C   ,  along  with  several  other  quantities,  must  be  computed 

AA 

at  station  K  in  order  to  get  A2T. .   Now,  Cn   is  defined  by 

3      K  XX  * 

equation  (3.4a),  and  is  computed  from 


„n  1   "    i  ■,    c   i+n    c   |i-n| 

cxx  =  4  I    w'c  ^    w'?  +  u  w,?   ) 


where 


nc  = 


1  if  n  >  0 
0  if  n  =  0 


"See  Appendix  B  of  Reference  [30 
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pc  =  f  if  i  ¥  n 

\2  if  i  =  n 
Since  C    must  be  computed  at  node  K,  w,   must  be  known 

X..X.  s 

at  whole  station  K  for  each  value  of  n.   As  evidenced  by 
(3.18b),  the  first  derivative  of  w   can  be  expressed  only 
at  the  half  stations.   Consequently,  an  interpolation  of 
some  form  must  be  utilized  in  order  to  define  w,j}  at  a 
whole   station  in  terms  of  w   at  nearby  whole  stations. 
The  same  situation  exists  for  several  of  the  other  non- 
linear terms.   The  interpolations  used  to  define  the  other 
derivatives  falling  in  this  category  are  developed  in 
Appendix  B.   One  example  of  these  interpolations  is  the 
expression  for  w1?,.  at  whole  station  K 

(W'  g  }K  =  lid  U0  (WK+1  "  WK-1>  "  (wK+2  "  Wk-2H 


Location  of  Rings 

In  the  interest  of  simplicity,  the  assumption  was 
made  that  rings  could  be  located  at  whole  stations,  but  not 
at  half  stations.   As  an  aid  in  expressing  the  ring- 
associated  terms ,  let  the  set  H  and  the  symbol  6     be 

KHj 

defined  as  follows: 


H  =  {H.  such  that  H.  =  K,  where  K  is  the  whole 

3  D 

station  location  of  the  j    ring;  j  =  1,2,...,N  } 

(3.19a) 
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1  if  K  =  H, 
6„„  =  4  :  (3.19b) 

KHj    \o    if  K  ?   H. 


IV.   DIFFERENCE-DIFFERENTIAL  EQUATIONS 

Equations  of  Motion 

In  this  section  the  difference-differential  equations 
which  constitute  the  numerical  model  are  obtained  by  re- 
placing  the  ^-derivatives  in  the  governing  equations  by 
the  modified  differences  of  (3.18). 

When  the  response  is  not  symmetric  with  respect  to  the 
midspan,  the  finite  difference  grid  is  arranged  so  that  the 

whole  stations  K  =  1  and  K  =  M  are  located  at  £  =  0  and 

2 
E,   =  i  ,  respectively.    The  total  number  of  whole  stations 

is  thus  equal  to  M.   The  half  stations  k  =  0  and  k  =  m  are 
imaginary  points  located  a  half  mesh  space  from  the  ends  of 
the  shell.   The  whole  stations  K  =  0  and  K  =  M+l  are  imag- 
inary points  located  a  full  mesh  space  from  the  ends  of 
the  shell.   Thus  there  are  M  whole  stations  and  m-1  half 
stations  on  the  shell. 

When  the  motion  of  the  shell  is  symmetric  with  respect 
to  the  midspan,  the  finite  difference  mesh  is  arranged  so 
that  there  is  a  whole  station  at  £  =  0  and  a  half  station 
at  £  =  i/2.   The  whole  station  at  £  =  0  is  assigned  the 
index  K  =  1,  the  half  station  at  £  =  i/2  is  assigned  the 
index  k  =  m,  and  m  =  M  is  now  the  number  of  half  stations 


2See  Figure  3.2. 
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as  well  as  the  number  of  whole  stations  located  on  the 
half  of  the  shell  between  £  =  0  and  £  =  \/2.      Thus,  the 
whole  station  K  =  M  is  located  at  £  =  \\    -  ^d . 

By  expressing  (3.6)  at  the  half  stations,  and  (3.7) 
and  (3.10)  at  whole  stations,  the  difference-differential 
equations  that  govern  the  response  can  be  obtained  in  the 

form 

N 

u   ,  +  dn(±_^)  vvxl  -[2  +  d  n   (— -)  -   I  d  n  (C.n 
k+1        2     K+l  2      j.  =  -^  4 

+  C5»  6KHj]  <   -   dn^    Vl   +  Vl  "  "  d'X1k     (3-20) 

Nr 

-  ^  *»  cc4  («i2  -  i)  -c5]  «KH.  t-»+1-w»> 

K=k=l, . . . ,M-1  (or  M) ;  n=0,l,... 


3 
The  rings  were  assumed  to  be  located  only  at  whole 

stations.   However,  in  equation  (3.20),  the  term 

N 

?  2      7 
-    Z  d  n   (C  n   +  C  )  6 
j  =  l        4       5    KHj 

is  multiplied  by  u,  ,  a  displacement  which  is  defined  at  a 
half  station  located  one-half  of  a  mesh  space  from  the  ring 
at  station  K  =  H • .   Since  the  axial  displacements  are  gener- 
ally small  compared  to  the  radial  displacements,  this 
inconsistency  is  probably  not  significant,  especially  if  the 
grid  points  are  close  together.   Alternatively,  the  above 
term  could  be  replaced  by 

Nr 

-  Z   *4d2n2(C  n2  +  C  )(6     +  6       ) 
-j--^         4       5    KH-     K+l  ti. 

Then  half  of  the  effect  of  a  ring  is  felt  at  station  k  if 
there  is  a  ring  at  station  K  or  K+l. 


77 


N 

,l~v.       n  #l±v%      n         rn  ,2    2  „r     2    2  ,_    x       on 

(— )    vK+1   -   dn    (— )    uk   -     [1-v+d   n      -      Z   d   n    (C3<SKH,)-I   VK 

+  dn    (    2    )    uk_x   +   —     v^   =   -   d    X2K 

Nr  «  (3.21) 

-I      d2nc3    (en-1)    6RH     wR 

3=1  D 

K=k=l,. . . ,M-1     (or   M) ;    n=0,l,.. 

dl^K   =    60n    PK   +   d2^K    +   d3WK+2    +   Vk+1    +   d5    WK 

,      n  n  n  n  n 

4    K-l  3    K    l  6  r  ^ 


K  <=K  ^UK 

N. 


n  r  n 


+    n         +      Z      6  R.,  (3.22) 

6K         j=l      KHj       D 

K=l , . . . ,M;    n=0 ,1 , . . . 


Here  d, , . . . ,d   are  known  quantities  which  depend  on  n,  c, 

jjmu  -v-.ni-,nn  n  n  ,      n 

a,    and   d.      The   terms    Al    ,    X2    ,    n      ,    i\      /    nr0    ,    and    n        are 

k    k   ek   eK   c«K      ^eK 

calculated  as  explained  in  Appendix  B.   In  (3.20)  and 
(3.21),  the  upper  limit  M  of  the  values  assumed  by  K  applies 
only  when  the  response  is  symmetric  with  respect  to  the 
midspan.   The  ring  term  R . ,  obtained  from  (3.11)  and  (3.18), 
is  given  by 

Rj  =  d6WK+l  +  d7WK  +  d6WK-l  +  d8(Uk  "  Uk-1>  +  dgVK 

K=H  ;  j=l,2, . . . ,N        (3.23) 
J  r 

Here  d  , ...,d   are  also  known  quantities  that  depend  on  n, 
6       9 

d,  and  the  ring  coefficients  C  , . . . ,C_ ,  and  e.   The  expres- 

1       5 

sions  for  d. , . . . ,d_  are  given  in  Appendix   C. 
i      y 
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Boundary  Conditions 

Casting  the  boundary  conditions  in  finite  difference 
form  is  conceptually  simple,  but  the  process  involves 
development  of  a  few  interpolations.   The  interested 
reader  will  find  the  development  of  the  interpolations  in 
Appendix  B. 

Simply  Supported.   The  boundary  conditions  (3.12) 
allow  an  axisymmetric  rigid-body  translation  along  the 
^-coordinate .   Consequently,  when  only  the  conditions  of 
(3.12)  are  observed,  the  numerical  solution  is  not  unique. 
In  order  to  obtain  a  unique  solution,  the  axial  displace- 
ment of  the  axisymmetric  mode  must  be  set  to  zero  at  some 
arbitrary  axial  station.   To  facilitate  a  numerical 
solution,  the  node  at  £  =  i  where  K  =  M  was  chosen.   Since 
the  axial  displacement  is  not  defined  at  the  whole  stations, 
interpolation  is  required  to  express  this  condition.   The 
interpolation  is   developed  in  Appendix  B.   Substituting 
the  appropriate  equations  from  (3.18)  into  (3.12)  at  K  =  1 
yields 

w£  =  0  (3.24a) 

w§  -  2w^  +  Wq  =  0  (3.24b) 

vn  =  0  (3.24c) 

u5  "  US  +  d<CxXl  +  ^XX^  =  °  (3'24d) 
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and   at   K   =   M  yields 


wJJ  =    0  (3.25a) 

WM+1    "    2wM   +   WM-1    =    °  (3.25b) 

vJJ  =    0  (3.25c) 

um   "    u^    +   d(C^        +   C5XX    )  =    0                                          (3.25d) 

M  M 


In  order  to  tie  down  the  axial  displacement  of  the  axi- 
symmetric  mode. the  following  expression  is  used  at  the  end 

i-  >4 

u°  +  2u°    -  -  u°    =  0  (3.26) 

m     m-1    3   m~2 

where  the  half  station  k  =  m  is  located  at  E,   =    i  -   \& . 

Clamped.   The  finite  difference  formulation  of  boun- 
dary conditions  (3.13b)  and  (3.13c)  is  straightforward 
requiring  only  the  substitution  of  the  appropriate  expres- 
sions from  (3.18).   However,  the  discretization  of  (3.13a) 

and  (3.13d)  requires  the  utilization  of  interpolation  form- 

4 
ulas.    The  difference  expressions  for  the  clamped  boundary 

conditions  are 

Uq  +  2u^  -  j  u^  =  0                             (3.27a) 

v^  =  0  (3.27b) 

w^   =  0  (3.27c) 


See  Appendix  B. 
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Wq    -    9v>2    +  W3   =    0  at   K   =    1  (3.27d) 


and 


^  +  2Ci  -  k-2  -  °  (3-28a> 

v£  =    0  (3.28b) 

M 

wJJ  =    0  (3.28c) 

w[!   0    -    9w£   .    +  w",    =0  at   K  =  M  (3.2  8d) 

M-2  M-l  M+l 

Response  Symmetric  to  the  Midspan.  Substitution  of 
the  appropriate  expressions  from  (3.18)  into  (3.14)  leads 
to 

u^  =  0  (3.29a) 

m 


VM  "  v£+l  =  °  (3.29b) 

WM  "  WM+1  =  °  (3'29C) 


WM-1  "  "M+2  =  °  (3-29d) 
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CHAPTER  IV 

SOLUTION  OF  THE  NUMERICAL  MODEL 

I.   SOLUTION  PROCEDURE 

There  are  various  methods  that  could  be  used  to  solve 
the  differential-difference  equations  that  comprise  the 
numerical  model.   Several  of  them  are  examined  in  this 
chapter.   The  algorithm  that  appears  to  be  the  most  promis- 
ing is  a  combination  of  the  Gauss  elimination  method  and 


the  iterative  Newmark  beta-method.   The  solution  procedure 

n 


is  carried  out  in  the  following  manner.   Suppose  that  w„, 


n      . n 
quantities  wK  and  wv  at  t0  are  then  calculated  for  all  K 

and  n  using  the  Newmark  beta-method.   Knowing  wK  at  x0,  the 

terms  Alv  and  A2„,  which  depend  only  on  wv ,    are  computed; 


w__,  wT. ,  u,  ,  and  vv   are  known  at  time  t,  for  all  values 

K     J\     K  ^  -L 

of  K,  k,  and  n.   The  solution  is  sought  at  time  t  ,  a  small 

increment  later  than  t, .   In  the  first  iteration  for  the 

solution  at  x_.  w   at  xn  is  taken  equal  to  w£  at  t, .   The 
2    K      ^ 

n      .n 

K  and  wK  at  T- 

TK  au  l  o 

and  u,  and  v„   are  obtained  by  solving  equations  (3.20)  and 
(3.21)  with  the  Gauss  elimination  method.   Next,  having 
computed  estimates  for  w„,  wK,  u,,  and  vK  at  time  t  ,  the 
acceleration  wK  at  t   can  be  calculated  from  equation  (3.22) 
for  all  K  and  n.   The  accelerations  of  the  two  consecutive 
iterations  are  compared  to  check  for  convergence.   The 
iteration  is  continued  until  convergence  is  achieved. 
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This  procedure  is  repeated  for  many  time  steps  to  compute 
the  displacements  as  a  function  of  time. 

II.   SOLUTION  OF  THE  AXIAL  AND  CIRCUMFERENTIAL  EQUATIONS 

Applicable  methods 

The  axial  and  circumferential  ordinary  differential 

equations,  (3.6)  and  (3.7),  constitute  a  linear  boundary 

n         n  n 

value  problem  in  the  two  unknowns  u  (£)  and  v  (£)  when  Al  , 

A2  ,  and  w   are  known.   In  the  development  of  the  numerical 

model,  these  equations  were  replaced  by  a  system  of  simul- 

n 
taneous  algebraic  equations  (3.20)  and  (3.21).   If  Al,, 

A2  ,  and  wK  are  known,  these  equations  can  be  solved  in- 
dependently of  the  radial  equations  (3.22).   Either  a 
direct  elimination  procedure  or  an  iterative  procedure  is 
normally  employed  to  solve  systems  of  algebraic  equations 
of  this  type  [35,  36].   An  example  of  the  direct  procedure 
is  the  Gauss  elimination  method.   This  method  has  been  used 
by  several  investigators  to  solve  static  shell  problems 
[23,  30,  37].   The  Gauss-Siedel  method  is  a  popular  itera- 
tive procedure. 

The  Gauss  elimination  method  and  the  Gauss-Siedel 
iteration  procedure  were  compared  in  an  attempt  to  determine 
the  best  of  the  two  methods  for  the  problem  under  considera- 
tion.  Computer  routines  using  both  methods  were  written 
to  solve  (3.20)  and  (3.21).   Solutions  to  some  simple 
linear  problems  were  generated  by  both  routines,  and  the 
computed  solutions  were  compared  with  the  exact  solutions 
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to  the  differential  equations.   The  Gauss  elimination  pro- 
cedure required  significantly  less  time  on  the  computer 
than  did  the  iterative  procedure  for  the  same  degree  of 
accuracy.   When  both  methods  used  the  same  execution  time, 
the  elimination  method  provided  a  considerably  more  accurate 
solution.   Consequently  the  Gauss  elimination  method  was 
selected  for  the  algorithm. 

Gauss  Elimination 

To  facilitate  the  development  of  the  Gauss  elimination 
procedure,  the  axial  and  circumferential  equations  are 
converted  to  matrix  equations.   First,  let  Al,  and  A2K  be 
defined  by 

Nr 
Alk  =  d'  Xlk  +  -\     6KH.  dn2[C  (en2-l)  -  Cg] <w£+1  -  w£) 

(4.1a) 

Nr 

A2^  =  d2  A2^  +   E   6TZU   d2n  C^(en2-1)  w£  (4.1b) 


Equations  (3.20)  and  (3.21)  can  now  be  written  in  the  form 

Anzn    +  Bnzn  +  Cnzn    =  a11  (4  2) 

A  ZK+1    bKZK  +  L  ZK-1    gK  [  ' 


where  n  =  0,1,...;  and  K  =  1,2,...,M-1  or  M,  depending  on 
whether  the  response  is  unsymmetric  or  symmetric  with 
respect  to  the  midspan.   For  the  unsymmetric  response,  M  is 
the  number  of  whole  stations  between  £  =  0  and  E,   =    \  .      If 
the  response  is  symmetric  with  respect  to  the  midspan,  M 
is  the  number  of  whole  stations  located  between  £  -  0  and 


84 


£=i/2.   The  matrices  An,  b£,  and  C   are  2x2  matrices  whose 
elements  are  given  in  Appendix  C,  and  zK  and  gK  are  defined 
by 


n 
ZK  = 


u. 


n 


v 


n 
K 


(4.3) 


n 


-Al 


n 


~^K 


(4.4) 


For  a  fixed  value  of  n,  (4.2)  represents  a  system  of  2M-2 
(or  2M)  equations  in  2M+2  (or  2M+4)  unknowns  for  unsymmetric 
(or  symmetric)  response  as  shown  in  Figure  4.1.   The  bound- 
ary  (or  symmetry)  conditions  provide  the  remaining  neces- 
sary equations 

Recursion  Relations.   The  Gauss  elimination  recursion 
relation  for  equations  having  the  form  of  (4.2)  is  [30,  37] 


n      n   n       n 
ZK  '  ~FK  ZK+1  +  XK 


(4.5a) 


where  P„  is  a  2x2  square  matrix,  and  xv  is  a  2x1  vector. 


'K 


n 


Thus,  the  recursion  relation  for  zK_-,  is 


n        n    n    n 
ZK-1      FK-1  ZK    XK-1 


(4.5b) 


Substitution  of  (4.5b)  into  (4.2)  results  in 

n     ,n    n  n    -1  n   n     rn    n  n    -1   n   n   n 
ZK  =  "[BK  "  C  PK-ll   A   ZK+1  +[BK  "  C  PK-J    ^K"C   XK-13 


85 


o 

o--x — o — 


2 

-o X- 


(AXIAL    EQUILIBRIUM 
1       SATISFIED  AT 
3^*\HALF  STATION  NUMBE 
>     x   ^  {—%. — o — x— o     *      o-( 


^CIRCUMFERENTIAL 
\    EQUILIBRIUM 
SATISFIED    AT 


'  \    4     \   r^\WHOLE  STATION  NUMBER \      M_l  A 

V     V     V      V  /       t^ 


INDEX      K    OF 
EQUATION    (4.2) 

(a  )UNSYMMETRIC     RESPONSE 


M 


M-l 


C  =  o 


I  2 

O-  -K r> X^-O x— ( 

0 


S>    £v> 


*  *  /2       f  AXIAL    EQUILIBRIUM 
m-l        m      ^J      SATISFIED   AT 
f_  -*-<  HALF  STATION  NUMBER 

N— o^x — o— X-— o xV. 

M  +  |    /^CIRCUMFERENTIAL 
EQUILIBRIUM 
SATISFIED    AT 


M-l      M 


WHOLE    STATION    NUMBER 


M-*-INDEX    K    OF    EQUATION    (4.2) 


M- 


Ib )  SYMMETRIC     RESPONSE 


FIGURE 


4.1 


RELATIONSHIP    BETWEEN  MESH  NODES 

AND 
INDEX    K    OF  EQUATION  (4.2) 


86 


Comparing  the  above  equation  with  (4.5a)  reveals  that 

p£  =  [Bn  -  cV   ]_1  An  (4.6) 

K      K       K-l 


and 


*K=  [BK  "  ^K-ll"1^  "  ^  4-S  <*•?> 


In  solving  (4.2),  the  following  steps  are  performed 
(1)  the  P^  matrices  are  computed  by  recursion  relation 
(4.6);  (2)  the  x£  vectors  are  calculated  from  (4.7);  and 
(3)  the  solution  vectors  z£  are  computed  in  a  backward 
sweep  by  recursion  relation  (4.5a).   Each  of  the  steps 
will  now  be  examined  in  more  detail. 


Computation  of  PK  Matrices.   The  general  matrix  PK 


depends  on  P^_i.   Hence,  P^  must  be  determined  first.   This 
matrix  is  determined  from  application  of  the  boundary  con- 
ditions.  The  boundary  conditions  which  uniquely  determine 

n 
the  Pj  matrix  are  given  by  equations  (3.24c)  and  (3.24d) 

when  the  shell  is  simply  supported  and  by  (3.27a)  and 

(3.27b)  when  the  shell  is  clamped.   From  the  two  appropriate 

boundary  equations  and  the  axial  equilibrium  equation   at 

K  =  1  an  equation  is  obtained  of  the  form 

n    ,-^n   n     n  .  .  _ . 

Zl  =   l  Z2  +  xl  (4,8) 

The  elements  of  the  matrix  P,  for  both  the  simply  supported 
and  clamped  shell  are  listed  in  Appendix  C.   By  using  the 
recursion  formula  (4.6),  the  remaining  PK  matrices  are 
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computed  for  all  of  the  participating  modes  and  for 
K  =  2,  3,  ...,M-2.   When  the  response  is  symmetric  with 
respect  to  the  midspan,  K  =  2,  3,  ...,M-1.   These  P 
matrices  do  not  change  with  either  the  load  or  the 
response . 


n 


Computation  of  x„  Vectors.   When  the  shell  is  simply 

1 
supported  the  vector  x?  of  (4.8)  is  given  by 


n 

xi = 


Al?  +  *<(£    +  C?    ) 


2  2  1-v 
/[l  +  d  n  (— -)] 


(4.9) 


and  when  the  shell  is  clamped,  by 


n 
x  = 

1 


Al. 


n 


/[4  +  d2n  i1-^)) 


(4.10) 


Since  in  both  cases  x,  contains  nonlinear  terms  in  the 
radial  deflection,  x-i  is  not  constant,  but  instead,  varies 
with  the  load.   The  vector  gK,  defined  by  (4.4),  also 
changes  with  the  load.   Thus,  the  vectors  xK  and  gK  must 

be  recomputed  many  times  as  a  solution  is  generated;  Xi 

n 
from  (4.9)  or  (4.10),  and  the  other  xK  vectors  from  (4.7). 


1  .  .        n  .    n 

The  fact  that  the  boundary  condition  on  v„  is  V]_  =  0 

makes  it  unnecessary  to  introduce  vQ  into  the  solution  and 

to  consider  the  circumferential  equilibrium  equation  at 

point  K  =  1 . 
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Back  Substitution.   The  solution  vectors  zR  are  ob- 
tained through  back  substitution,  that  is  by  sequential 
solution  of  equation  (4.5a).   However,  relation  (4.5a) 
requires  a  "starting"  z-vector. 

The  boundary  conditions  at  £  =  i  provide  the  infor- 
mation necessary  to  compute  this  vector  when  the  response 
is  unsymmetric  with  respect  to  the  midspan  of  the  cylinder. 
The  boundary  conditions  are  given  by  (3.25c)  and  (3.25d) 
when  the  shell  is  simply  supported  and  n  ^  0 .   When  the 
shell  is  clamped,  or  when  the  shell  is  simply  supported 

and  n  =  0,  the  corresponding  boundary  conditions  are  ex- 

2 
pressed  by  (3.28a)  and  (3.28b).    By  substituting  the 

boundary  conditions  into  (4.2)  with  K  =  M-l  and  by  making 

use  of  (4.5a)  at  K  =  M-2,  the  following  result  is  obtained 

for  the  solution  at  K  =  M-l. 

n      rTTn    ^n_n   ,  ,,.n    -.n   n   ,  .  .    n  ..  . 

ZM-1  =  [H   "  °  PM-2][f   '  °   XM-2]  (4'1:L) 

where  H  ,  0  ,  and  f   are  given  in  Appendix  C  for  the  simply 

supported  and  clamped  boundary  conditions. 

When  the  response  is  symmetric  with  respect  to  the 

midspan,  the  starting  vector  zM  is  determined  by  using  the 

symmetry  conditions  given  by  (3.29a)  and  (3.29b).   If 

(3.29a)  and  (3.29b)  are  substituted  into  the  circumferential 

equilibrium  equation  at  K  =  M  and  the  result  combined  with 

(4.5a)  at  K  -  M-l,3  the  result  is 


2 
Equations  (3.26)  and  (3.28a)  are  identical  when  n=0. 

3 
The  axial  equilibrium  equation  is  identically  satis- 
fied at  the  midspan. 
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1         n 

zJJ  =  [I  +  Rn  Pm_;l]~  [^  *M-1  +  sn]  (4.12) 

n       n       , 
where  I  is  the  identity  matrix  and  R  and  s   are  given  in 

Appendix  C. 

Once  equations  (4.11)  or  (4.12)  have  been  evaluated 

for  each  n,  the  remaining  z-vectors  are  computed  in  sequence 

using  (4 .5a) . 

III.   TIME  INTEGRATION 

Four  basic  methods  were  considered  for  carrying  out 
the  time  integration  of  the  radial  equations  of  motion. 
They  were  the  Runge-Kutta  methods,  the  explicit  methods, 
the  implicit  methods,  and  the  Newmark  beta-method. 

The  Runge-Kutta  method  [38]  has  been  used  by  several 
investigators  to  compute  the  nonlinear  dynamic  response 
of  cylindrical  shells  [11,  20,  22].   However,  they  included 
only  a  few  modes  in  the  response  and  did  not  employ  finite 
differences  along  the  axial  coordinate.   Consequently, 
their  systems  of  equations  were  much  smaller  than  the  system 
to  be  solved  here.   For  the  large  system  of  equations  of  the 
present  analysis,  the  Runge-Kutta  method  would  require 
prohibitively  large  amounts  of  computer  storage  space  and 
excessive  computing  time. 

The  explicit  methods  [39,  40]  all  have  the  common 
feature  that  the  unknowns  can  be  calculated  explicitly 
without  solving  a  large  system  of  simultaneous  equations 
each  time  the  solution  is  advanced  one  time  increment. 
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Nonlinear  equations  can  be  solved  without  iteration  since 
the  nonlinear  terms  contain  only  products  of  previously 
computed  quantities.   However,  the  explicit  methods  are 
numerically  unstable  if  the  time  increment  is  too  large. 
While  the  largest  permissible  time  increment  can  be  pre- 
dicted theoretically  for  linear  equations,  the  permissible 
time  increment  for  nonlinear  equations  has  not  been  estab- 
lished as  yet.   Since  a  numerical  instability  might  easily 
be  mistaken  for  a  physically  unstable  response,  the  ex- 
plicit methods  were  ruled  out. 

The  implicit  methods  [39,  40,  41]  all  utilize  time 
averaged  finite  difference  approximations  of  the  spatial 
derivatives.   As  a  result,  the  difference  equation  at  any 
spatial  location  contains  several  unknowns.   The  advance- 
ment of  the  solution  by  one  time  increment  is  accomplished 
by  solving  a  system  of  simultaneous  algebraic  equations. 
If  the  differential  equations  are  nonlinear  in  the  spatial 
derivatives  of  the  dependent  variables,  as  in  the  case  for 
the  equations  under  study,  the  simultaneous  algebraic 
equations  are  usually  nonlinear.   These  nonlinear  equations 
must  be  solved  by  an  iterative  procedure  that  solves  a 
sequence  of  linear  equations.   The  implicit  methods  have 
the  advantage  of  being  numerically  stable  when  applied  to 
linear  problems.   However,  since  an  iterative  procedure 
must  be  used  in  the  implicit  methods  to  solve  the  nonlinear 
equations,  the  decision  was  made  to  use  an  iterative  inte- 
gration procedure. 


91 


The  iterative  integration  procedure  selected  as  the 
most  suitable  for  carrying  out  the  time  integration  is  the 
Newmark  beta-method.   References  [42]  and  [43]  contain  a 
detailed  description  of  the  method  and  a  theoretical  deter- 
mination of  the  stability  and  convergence  limits.   A  typical 
application  of  the  method,  as  well  as  references  to  several 
other  investigations  that  employed  the  method  may  be  found 
in  Reference  [44]. 

Integration  Procedure 

The  numerical  integration  procedure  begins  with  the 

assumption  that  the  displacements,  u^,    vK,  and  w  ;  the 

•  n  ..n 

velocities,  wK;  and  the  accelerations,  wK,  are  known  at 

time  t  .   The  corresponding  physical  quantities  may  be 
computed  at  time  x_  =  x.  +  At,  where  At  is  the  time  incre- 
ment, by  repetition  of  the  following  steps: 

1.  Assume  that  wk(t  )  =  wk(t  )  for  the  first 

iteration.   For  subsequent  iterations,  assume 
that  wv(t9)  =  w1k(t  )  where  wlK  is  computed 
in  step  6 . 

2.  Calculate  the  radial  displacements  and  velocities 

at  time   t   by  the  relations 

^K(T2)  =  *K(T1)  +  1^t^K(Tl)  +^K(T2)]  AT         (4.13) 

WK(T2)  =  wK(TiJ  +  ^k(tx)At  +  ft-S)  w^(t1)(At)2  (4.14) 

..n         2 
+  3wk(t2) (At) 
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The  significance  of  the  parameter  8  will  be 

discussed  in  a  later  paragraph. 

n 


3.   Evaluate  the  nonlinear  terms  Alv  and  A2V  from 


(4.1)  . 

4.  Compute  uk  and  vK  by  Gauss  elimination. 

r    t,      i        n    n     n        n 

5.  Evaluate  nn  ,  rir  ,    nro  t    and  ru 

6K    4    ^K        9K 

6.  Compute  the  radial  accelerations  from  equation 

(3.22).   Let  wl„(Tn)  denote  the  values  computed 
here  to  distinguish  them  from  wk(t  )  of  step  L 

7.  Compare  w1^(t?)  with  wK(x2).   If  their  relative 

difference  is  larger  than  a  prescribed  value  z, 
the  convergence  criterion,  repeat  steps  1-7. 
If  the  difference  is  smaller  than  e,  the  solu- 
tion is  said  to  have  converged. 
This  procedure  can  be  repeated  as  often  as  desired 
in  order  to  determine  the  response  of  the  shell  as  a  func- 
tion of  time. 

The  Parameter  Beta.   Equations  (4.13)  and  (4.14)  are 
the  quadrature  relations  of  the  Newmark  beta-method. 
Newmark  [4  3]  has  shown  that  the  choice  of  a  particular 
value  of  6   is  related  to  the  shape  of  the  acceleration- 
time  curve  during  the  time  increment  At .   Acceleration- 
time  curves  associated  with  6  =  1/4,  1/6,  and  1/8  are  shown 
in  Figure  4.2.   For  this  study,  8  was  taken  equal  to  1/6 
corresponding  to  a  linear  variation  of  the  acceleration 
during  the  time  increment. 
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wl 


6=1/6 


Figure  4.2 
ACCELERATION  VS.  TIME  FOR  6  =  1/4,  1/6,  AND  1/8 

The  stability  and  convergence  limits  of  the  Newmark 
beta-method  are  given  by  critical  values  of  the  ratio  of 
the  time  increment  At  to  the  shortest  natural  period  T  of 
the  finite  difference  model.  If  the  ratio  At/T  is  larger 
than  the  stability  or  convergence  limit,  the  procedure  is 
unstable  or  does  not  converge.   The  stability  limits, 
derived  by  Newmark  [43]  for  a  linear  oscillator,  are  given 
in  Table  I  for  several  values  of  3.   When  $  is  larger  than 
one-eighth,  the  convergence  limit  is  more  stringent  than 
the  stability  limit. 

The  problem  under  study  is  governed  by  a  set  of  non- 
linear equations,  and  the  shortest  natural  period  is  not 
easily  obtained  since  the  natural  frequencies  of  a  non- 
linear system  vary  with  the  amplitude  of  the  motion. 
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Nevertheless,  the  shortest  natural  frequency  of  an  un- 
stiffened  shell,  as  predicted  by  the  linear  theory,  will 
serve  as  a  useful  guide. 

TABLE  I 
STABILITY  AND  CONVERGENCE  LIMITS 

3  0       1/12     1/8     1/6      1/4 

Stability  Limit, 

At/T  =  0.318     0.389   0.450   0.551   infinite 

Convergence  Limit, 

At/T  =  infinite   0.551   0.450   0.389    0.318 


For  the  Fourier  modes  of  interest,  the  linear  equa- 
tions reveal  that  the  axisymmetric  mode  (n=0)  has  the 

4 
shortest  period.    If  the  shell  is  simply  supported, 

nr  =  n   =  0  and  n   =  -(1-v  )w  .   Hence,  from  (3.10)  ,  the 

linear  differential  equation  for  the  axisymmetric  mode  is 

seen  to  be 

12   0  2    0     .0 

12  a      ™'EEEE+    ^1_V  )  w   =  "w  (4.15) 

Approximation  of  the  spatial  derivatives  by  central  differ- 
ences leads  to 

lrv2/0  ,0  0  ,0  0.  ,  .    .-, 

TJ  \     <wK+2    "    4wK+l    +   WK    "    4WK-1    +   WK-2>  (4*16) 

2,        0 
+  (1-v^)    wK   =    -wR 


4In  Chapter  V,  it  will  be  shown  that  this  is  true. 
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Assume  a  solution  to  (4.15)  of  the  form 

w°  =  sin  (ftQ  t)  sin  (32i)  (4.17) 

where  ft    is  the  natural  frequency  of  the  axisymmetric 
mode  with  wavelength  2—  ,   if  the  shell  is  divided  into 
M-l  intervals  of  length  d,  it  follows  that 


?  =  (K-l)  d 
i    =  (M-l)  d 

Thus  equation  (4.17)  may  be  written 

w°  =  sin  (ft.  t)  sin  [qTT(K"1)  ]  (4.18) 

Oq  M_]_ 

The  natural  frequency  ft    is  now  the  frequency  associated 

Oq 

with  the  difference-differential  equations.   Substituting 
(4.18)  into  (4.16)  and  performing  several  algebraic  and 
trigonometric  manipulations  yield 

Q,n      =  —  —  [2  cos  (—2—)  -  8  cos  (rr-r)  +  6]  +  (l-v  ) 
Oq    12  d4  M-l  M-l' 

9 

The  maximum  value  of  ft~   occurs  when  -3—  =  1.   Hence,  the 

Oq  M_! 

highest  frequency  is  given  by 

ft^      =  i  "   +  l-v2  (4.19) 

0  M-l    3  d4 

Consequently,  the  shortest  natural  period  of  the  finite 
difference  model  is 


,2       o  —2 


-h 
T  =  2tt  (i  2L.  +  l-v2)  (4.20 

J  d 
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With  T  given  by  (4.20),  At  for  convergence  and 
stability  is  given  in  Table  I  as  a  function  of  3-   Thus 
taking  3  equal  to  1/6 , 

At  <  (0.389)  (2tt)  (t  \  +  J.-V  )  (4.21) 

J  d4 

insures  both  convergence  and  stability  of  the  algorithm 
in  computing  a  linearized  solution.   The  assumption  is 
made  that  the  same  reasoning  can  be  applied  to  the  non- 
linear equations;  that  is,  if  the  calculation  converges, 
it  is  also  numerically  stable. 

Convergence  Rate  and  Computational  Errors.   The  rate 
of  convergence  used  by  Newmark  is  defined  as  the  error  in 
the  derived  acceleration  divided  by  the  error  in  the  assumed 
acceleration.   Newmark  has  shown  that  for  a  linear  oscilla- 
tor, the  rate  of  convergence,  as  well  as  the  error  in  the 
amplitude  and  frequency  of  the  computed  motion,  depends 
on  the  size  of  the  time  increment.   These  relationships 
are  shown  in  Figure  4.3  for  3  equal  to  1/6.   Inspection 
of  Figure  4.3  reveals  that  in  order  to  keep  the  errors 
small,  the  time  increment  must  be  much  shorter  than  the 
convergence  limit  permits.   For  example,  if  errors  in 
amplitude  and  frequency  larger  than  two  percent  can  not  be 
tolerated,  then  At/T  must  be  less  than  0.12,  whereas  the 
convergence  limit  requires  only  that  At/T  be  less  than 
0.389. 
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Newmark  [42]  has  suggested  that  when  3  is  equal  to 
1/6,  a  stable  integration  procedure  for  nonlinear  equations 
can  be  obtained  by  requiring  that  the  assumed  acceleration 
and  the  derived  acceleration  agree  to  within  one  percent, 
or  less,  by  the  end  of  the  fourth  iteration.   If  conver- 
gence is  not  obtained  by  the  fourth  iteration,  the  time 
increment  should  be  reduced.   On  the  other  hand,  if  the 
accelerations  converge  in  one  or  two  iterations,  the  time 
step  should  be  increased  to  ensure  optimum  computing 
efficiency. 

V.   COMPUTER  PROGRAM 

A  computer  program  was  written  for  the  IBM  360  Model 
67  computer  to  generate  a  solution  to  the  governing 
difference-differential  equations  developed  in  Chapter  II 
and  Chapter  III.   Appendix  D  contains  a  description  and 
complete  documentation  of  the  program,  including  instruc- 
tions for  the  preparation  of  input  data  cards  and  other 
information  necessary  for  the  operation  of  the  program. 
The  program  is  dimensioned  to  accommodate  twenty  equally 
spaced  axial  stations  (M  <  20)  and  fifteen  Fourier  modes. 
It  will  handle  up  to  three  ring  stiffeners  located  so  as 
to  divide  the  shell  into  segments  of  equal  length.   Both 
simple  and  clamped  supports  can  be  handled.   The  basic 
output  consists  of  the  displacements  and  in-plane  forces 
at  specified  time  steps.   The  axial  location,  time,  and 
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magnitude  of  the  maximum  circumferential  stress  at  6  =  0 
are  also  printed. 

The  program  adjusts  the  size  of  the  time  increment  so 
that  the  number  of  iterations  always  falls  within  a  range 
specified  by  the  user.   The  number  of  iterations  required 
to  obtain  a  converged  solution  depends  on  the  convergence 
criterion,  £.   Convergence  is  assumed  to  have  occurred  when 

,,.n/+1    ,..n  * 

I (wK)     -  (w  )         n  =  0,1,... 

£   >  _ 


1+1  I 


Max  ||  (wg)    |  ,  |  (w£) 
where  I  is  the  iteration  number 


K  =  0,  . . .  ,M 


Calculation  of  the  Circumferential  Stress  at  a  Typical 
Location 

A  structure  designed  for  repeated  use  may  experience 
low  cycle  fatigue  cracking  and  even  total  failure  if  the 
stresses  repeatedly  exceed  the  yield  point  stress  of  the 
material.   Therefore,  one  criterion  for  a  safe  load  is 
that  the  stresses  which  develop  during  the  response  must 
never  exceed  the  yield  stress.   For  the  problem  at  hand, 
only  radial  loads  are  considered;  and  the  load  is  carried 
principally  by  circumferential  stresses.   Thus,  the  circum- 
ferential stress  is  a  close  approximation  to  the  maximum 
principal  stress.   In  the  present  program  the  circumferen- 
tial  stresses  are  computed  at  6  =  0°. 


^At  a  fixed  axial  location,  the  circumferential  stress  is 
not  necessarily  largest  at  0  =  0°,  but  since  the  assumed 
imperfection  modes  are  in  phase  at  6  =  0°,  it  is  a  reason- 
able indication  of  the  state  of  maximum  stress  in  the  shell. 
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At  a  point  on  the  midsurface  of  the  shell,  the  cir- 
cumferential stress  cjg  is  given  by 

a   =  r  nA  (4.22) 

6    1-v2   9 

At  either  the  inner  or  the  outer  surface  the  magnitude  of 

the  circumferential  stress  is  larger  than  the  stress  at 

the  midsurface,  the  additional  stress  being  due  to  bending 

of  the  shell.   The  circumferential  surface  stress  due  to 

the  circumferential  curvature  change,  denoted  a0lJ,  is  given 


by 


1  Eh  d2W  ..  ... 

a6B  "  ±  J  T  712  (4'23) 

a  do 


The  influence  of  the  meridional  curvature  change  on  the 
circumferential  bending  stress  is  neglected.   The  magni- 
tude of  the  largest  surface  stress,  denoted  I  cj   I,  is 

6s 


given  by 


°esl  "  IV  +  I'bbI  (4-24) 


Substituting  (4.22)  and  (4.23)  into  (4.24)  leads  to 

Ia0sl  =.^%|  n0  |  *l*.<*,l£*\  (4.25) 

6S     1-v2     y      2    a2    d62 


On  substituting  (3.2c)  and  (3.3b)  into  (4.25),  the  result 
is 

_  00  t  00 

i    i      Ei^n       n  i    !    ir,2n 


£   nQ  cos  n6   +  —  Ea   Z   n  w   cos  n8  (4.26) 


0S     1-v2    n=0   6  2 
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When  6=0° 


CD 

a 


Z   n"  |  +  i  Ea|  Z   n2wn |     (4.27) 


6S  0=0°    1-v2 


Let  nfl  be  defined  by  the  relationship 

I CJ   |      =  — 5J  n  (4.28) 

6S  0=0°    1-v2   6 

Then  from  (4.27)  and  (4.28)6 

CO  CO 

nQ  =  I  Z   n?  I  +  J-  a(l-v2)  I  Z   n2wn|  (4.29) 

n=0  n=0 


At  a  given  axial  station,  the  value  of  nfi  is  thus  propor- 
tional to  the  stress  at  0  =  0  .  Equation  (4.29)  can  also 
be  expressed  in  finite  difference  form  as 


oo 

n 


CO  CO 

Z      n£    |    +   i     a |    Z      n2w£ |  (4.30) 


6K  n=0      6K  2        '   n=0 


The  maximum  circumferential  stress  at  8=0  ,  nfl    ,  occurs 

MAX 
at  some  grid  point  K.   Hence 


lor  M)l 


nQ     =  Max  •{  nQ  ,  K  =  2,3,...,M-1  (or  M)y       (4.31) 
9MAX        I  6K 


For  convenience,  we  will  refer  to  this  as  the  maximum 
circumferential  stress. 


6  2 

In  the  computer  program,  the  quantity  1-v   in  (4.29) 

was  taken  equal  to  unity.   The  error  introduced  is  cer- 
tainly no  greater  than  the  errors  involved  with  the  other 
approximations.   If  the  program  is  revised  to  compute  the 
exact  stresses,  this  factor  should  be  reinserted. 
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Verification  of  Computer  Program 

As  indicated  in  Chapter  I,  there  apparently  are  no 
previous  solutions  for  the  nonlinear  dynamic  response  of 
a  cylindrical  shell  that  can  be  used  to  check  all  of  the 
features  of  the  present  computer  program.   Furthermore, 
the  program  cannot  be  readily  simplified  or  altered  to 
correspond  with  the  previous  analyses.   Consequently,  only 
a  linearized  version  of  the  program  was  verified  through 
comparison  with  other  solutions  known  to  be  correct.   It 
might  be  argued  that  some  of  the  new  features  of  the  pro- 
gram are  of  minor  significance,  and  that  the  results  of 
some  of  the  earlier  analyses  could  be  used  to  test  the  pro- 
gram.  However,  if  discrepancies  were  discovered,  it  would 
be  difficult  to  decide  whether  to  attribute  them  to  an 
error  in  the  program  or  to  differences  in  the  analyses. 

To  check  the  capability  of  the  program  to  compute  a 
fully  nonlinear,  nonaxisymmetric  solution,  a  typical 
problem  was  run  through  several  time  steps,  and  a  complete 
sequential  printout  of  the  program  variables  was  obtained. 
The  printout  was  then  checked  against  a  hand  calculated 
solution. 

To  date,  those  features  of  the  program  that  compute 
the  contributions  from  the  rings  have  been  checked  in  a 
cursory  manner.   One  computer  run  was  made  for  which  the 
shell  was  assumed  to  have  ring  stiffeners.   The  computed 
displacements  appeared  to  be  correct.   The  parts  of  the 
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program  that  compute  the  ring  contributions  are  isolated 
and  do  not  affect  results  obtained  for  unstiffened  shells. 
The  results  presented  in  the  next  chapter  also  serve 
to  verify  the  analysis  and  the  program  in  that  many  of  the 
results  agree  very  favorably  with  those  obtained  by  other 
investigators . 

Linear  Test  Cases 

Under  an  axisymmetric  load,  the  response  of  the  shell, 
as  governed  by  linearized  equations,  is  also  axisymmetric. 
The  correctness  of  the  program  in  computing  a  nonaxisymmetric 
response  was  checked  in  conjunction  with  the  verification 
of  the  nonlinear  features  of  the  program. 

Moving  Step  Load.   The  transient,  axisymmetric  re- 
sponse of  a  simply  supported  shell  subjected  to  a  moving 
step  pressure  discontinuity  has  been  determined  earlier 
by  P .  G.  Bhuta  [45]  using  a  series  solution  to  the  linear 
Donnell  equations.   A  subroutine  was  assembled  to  compute 
p(£,  t)  for  this  loading  condition,  and  a  numerical  solu- 
tion for  the  axisymmetric  response  was  computed  using  the 

program  with  all  nonlinear  terms  set  equal  to  zero. 

7 
Twenty-one  axial  stations  were  used,   and  convergence  to 

one  per  cent  in  a  maximum  of  four  and  a  minimum  of  three 

iterations  was  specified.   A  comparison  of  the  radial 


7 

The  present  program  is  dimensioned  to  handle  a 

maximum  of  twenty  axial  points. 
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displacements  at  the  quarter-span  station  is  shown  in 
Figure  4.4.   The  constants  for  the  example  are  also  given 
in  Figure  4.4.   The  computed  results  agree  fairly  well  with 
the  series  solution.   However,  by  the  end  of  the  calcula- 
tion a  detectable  error  has  developed. 

The  computed  displacements  are  given  in  Figure  4.4 
it  eviHry  tenth  time  step.   During  Lhe  time  that  the  load 
discontinuity  was  on  the  shell,  the  time  increment  remained 
small,  and  approximately  one  hundred  and  seventy  time 
stops  were  taken  to  compute  the  first  cycle  of  the  motion. 
Shortly  before  the  discontinuity  moved  off  the  shell,  the 
tune  increment  increased,  and  approximately  thirty  time 
steps  were  required  to  compute  the  second  cycle.   The 
ratio  At/T,  where  T  is  defined  by  equation  (4.20)  ,  reached 
a  maximum  value  of  .05.   From  Figure  4.3,  it  is  seen  that 
i_he  computational  errors  should  be  very  small  for  this 
value.   The  numerical  errors  are  apparently  accumulative 

nd  may  reach  signiiicant  magnitudes  in  lengthy  compu- 
tations.  In  subsequent  computations,  the  parameter  e  was 
set  equal  to  .001  to  reduce  the  accumulative  error. 

Uniform  Harmonic  Load.   The  axisymmetric  transient 
response  of  a  simply  supported  shell  subjected  to  a  loading 
cjiven  by 


p (£ ,x )  =  p   cos  Qi 


was  computed  and  compared  to  the  series  solution.   The 
response  predicted  by  the  linear  theory  is 
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FIGURE    4.4 
QUARTER-SPAN  REPONSE  TO  A  MOVING 
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31  4pQ 

w(£,x)  =  Z  Ty Ty—    (cos  fix  -  cos  Q    x)  sin  (2—3-) 

q=l,3,5,...  irq  (iT-iT)  q         ; 

p 
Thirty-one  axial  stations  were  used,   and  the  convergence 

parameter  was  set  equal  to  .001.   The  computed  solution 

at  every  10th  time  step  and  the  series  solution  are  shown 

in  Figure  4.5.   As  with  the  previous  test  case,  a  small 

error  developed  after  several  cycles  of  the  oscillation 

and  was  manifested  primarily  as  a  phase  difference  between 

the  computed  and  exact  responses. 

On  the  basis  of  the  above  tests,  it  was  concluded 

that  the  algorithm  provides  an  adequate  representation  of 

the  axisymmetric  transient  response  of  a  shell  subjected 

to  axisymmetric  loads. 


During  the  initial  phases  of  the  investigation,  the 
program  was  dimensioned  to  accomodate  31  axial  stations. 
The  dimensions  of  the  computer  program  were  later  reduced, 
and  the  program  listed  in  Appendix  D  will  handle  a  maximum 
of  twenty  axial  stations. 
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CHAPTER  V 

NUMERICAL  RESULTS 

I .   INTRODUCTION 

The  analysis,  the  numerical  model,  and  the  method 
of  solution  developed  in  the  preceding  chapters  are  appli- 
cable to  a  rather  broad  class  of  problems.   A  complete 
study  to  determine  the  significance  of  the  many  features 
incorporated  in  the  analysis  and  the  computer  program  is 
beyond  the  scope  of  the  present  work.   Hence,  a  problem 
that  has  been  considered  in  two  very  recent  investigations 
[19,  20]  was  selected  for  an  in-depth  study.   Treatment 
of  the  problem  by  the  algorithm  of  this  dissertation  ac- 
complished three  things:  (1)  increased  understanding  of 
the  solution  to  the  problem  was  developed;  (2)  the  feas- 
ibility and  practicality  of  the  new  approach  were 
established;  and  (3)  deficiencies  of  the  program  were  dis- 
closed. 

Statement  of  the  Problem 

The  problem  is  the  determination  of  the  dynamic 
response  of  an  unstiffened,  simply  supported   cylindrical 
shell  subjected  to  an  exponentially  decaying,  uniform, 
radial  pressure  pulse. 

Shell  Properties.   The  shell  was  assumed  to  be  made 
from  6061-T6  aluminum  and  to  have  the  following  properties 
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a    =  3 . 0  in 
L/a  =  2 

h/a  =  1/100 

-4         2.-4 
p    =  2.54  x  10   lb  -  sec   -  in 

E    =  107  psi 

v    =  0.3 


a    =42  ksi 

yp 

c    =  3.416  x  10"5 


where  a    is  the  yield  point  stress  of  the  material.   The 
yp         x  * 

nondimensional  damping  coefficient  c  corresponds  to  a 
value  determined  experimentally  by  Fung,  Sechler,  and 
Kaplan  [46].   Yao  [6]  used  the  same  coefficient  in  his 
dynamic  stability  study. 

Initial  Imperfection.   The  existence  of  a  small,  but 
finite,  initial  imperfection  in  each  mode  was  postulated. 
The  Fourier  coefficients  of  the  initial  imperfection  were 
assumed  to  be  of  the  form 

— n     .  n   .   ,  it  £ .  /  c  t  \ 

w  =  wi   sin  (— -)  (5.1) 

where  wi   is  the  amplitude  of  the  imperfection  in  mode  n. 
Except  for  one  special  case,  the  amplitudes  of  the  imper- 
fection were  all  taken  equal  to  l/1000th  of  the  shell 
thickness.   The  corresponding  nondimensional  amplitudes 
were  thus 

win  =  1.0  x  10~5 
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For  the  one  special  case,  the  imperfection  amplitudes  were 
taken  from  a  table  of  random  numbers. 

Load.   The  loading  was  assumed  to  be 


(5.2) 


where  P   is  the  peak  pressure  and  t   is  a  time  constant. 

o        c  c  o 

The  loading  is  expressed  nondimensionally  by 


( 

0                for    t    <    0 

p(x,t)  =\ 

__t_ 

P   e    to    for   t        0 
<  o                           — 

P(C,T)   = 


0       for  t  <  0 


x 
p  e   o  for  t  >  0 
*o  — 


(5.3) 


For  this  study,  t  >>1.   By  integrating  the  load  from 
t  =  -°°  to  t  =  +°°,  the  impulse  per  unit  area  is  found  to 
be 


I  =  PotQ  (5.4) 


II.   PREVIOUS  INVESTIGATIONS  OF  THE  PROBLEM 

Certain  aspects  of  the  problem  outlined  above  have 
been  the  topic  of  two  very  recent  articles.   Anderson  and 
Lindberg  [19]  investigated  dynamic  buckling  under  an 
exponentially  decaying,  uniform,  radial  pressure,  given  by 
(5.2).   Mclvor  and  Lovell  [20]  computed  the  dynamic  response 
for  a  uniform  radial  impulse,  a  limiting  case  of  the  loading 
considered  here  and  in  Reference  [19].   In  the  following 
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two  sections,  these  investigations  are  reviewed  in  some 
detail  in  order  to  lay  the  foundation  for  later  comparisons. 
The  simplified  equations  used  by  Anderson  and  Lindberg 
are  examined  in  considerable  detail  as  they  illustrate 
many  of  the  essential  features  of  the  problem. 

Anderson  and  Lindberg 's  Investigation 

Anderson  and  Lindberg ' s  analysis  identified  the  peak 
pressure  and  the  total  impulse  as  the  important  load  para- 
meters for  determining  whether  dynamic  buckling  will  occur. 
When  the  total  impulse  was  relatively  small,  they  found 
that  buckling  occurred  at  very  high  pressures;  in  fact,  so 
high  that  buckling  took  place  only  after  the  onset  of 
plastic  flow.   A  tangent  modulus  model  was  developed  for 
predicting  buckling  loads  in  this  regime.   For  loads  of 
large  total  impulse,  buckling  occurred  at  much  lower  pres- 
sures and  was  governed  by  an  elastic  model.   Both  of  these 
models  were  developed  for  the  prediction  of  dynamic  buckling 
loads  and  were  not  intended  to  predict  the  dynamic  response 
under  sub-critical  loads.   Their  elastic  model  is  of 
interest  here  and  will  be  used  as  a  standard  of  comparison. 

For  their  elastic  model,  they  assumed  the  motion  to 
consist  of  a  uniform  radial  displacement  plus  an  infinite 
number  of  nonaxisymmetric  flexural  modes.   The  flexural 
modes  were  required  to  satisfy  simply  supported  boundary 
conditons,  and  the  axial  variation  of  both  the  flexural 
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modes  and  the  initial  imperfections  was  taken  to  be  a  half 
sine  wave.   These  assumptions  were  satisfied  by 


0     „    n      ,tt£.       . 
w  =  w   +   Z   w   sin  (—2-)    cos  ni 

n=l 


for  the  radial  displacement  and 


—    ~    .  n   .   ,tt£ . 

w  =   Z   wi   sin  (— -)  cos  m 

n=  1 


for  the  initial  imperfection.   Finally,  the  coupling  between 
the  flexural  modes  plus  the  nonlinear  terms  in  the  equations 
of  the  axisymmetric  mode  were  neglected.   The  resulting 
equations  for  the  undamped  motion  were 

T 

w+w=pe°  (5.5) 

^o 


w  ,  1       2.    2         tt    .z,     (1-v    )  (tt/i)  2    o]    =   n   w 

+    N-T-  a    (n      +   — =■)    + \     —  -    n   w    J       .n 

.n  12  2  ~  2    ?  wi 

wi  i  .    2         tt    ,  ^ 

(n      +   -^ 

n=l,2,...  (5.6) 


Equation  (5.6)  can  be  written  in  the  form 

w  2  ,    n  O.w  20  ,0  /c-j\ 

+   n  (p         -   w    )    =nw  n=l,2,...  (5.7) 

.n  ^cr  . n 

wi  wi 


where 

pn  J_^a2    (n2   +  4^+a-v2)(VO^  (5.8) 

n  i  ,    2        it    .  z 

(n      +   —5-) 

1 
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The  nondimensional  static  pressure  which  will  cause  classi- 
cal buckling  of  the  n   mode  is  equal  to  p£   given  by  (5.8) 

+-  v» 
The  natural  frequency  of  the  n   mode,  Q    ,    is  obtained 

from  (5.7)  by  equating  w   to  zero,  from  which 

SI2   =   n2pn  (5.9) 

n     rcr 

Table  II  contains  the  natural  frequencies,  the  natural 

periods  T   =  2t\/Q.    ,  and  the  buckling  pressures  of  several 
r-        n       n 

of  the  modes  of  the  shell  under  consideration.   Note  that 
the  axisymmetric  mode  has  the  highest  natural  frequency, 
as  was  assumed  in  Chapter  IV.   Also,  mode  6  has  the  low- 
est buckling  pressure  and,  thus,  is  the  static  buckling 
mode . 


Hyperbolic  Response.   Anderson  and  Lindberg  [19] 
pointed  out  that  whenever  w  >p   ,  the  solution  of  (5.7) 
is  hyperbolic  in  nature,  and  the  amplitude  of  the  dis- 
placement grows  in  an  exponential  manner  with  time.   The 

closed  form  solution  of  equation  (5.5)  is 

2     _  T 
w  (t)  =  =-  [e   °  -  cos  t  +  —  sin  t]        (5.10) 

1  +  T  O 


when  w   and  w   are  initially  zero.   As  seen  from  Table  II, 
the  frequencies  of  modes  four  through  eleven  are  consider- 
ably less  than  ft-,  the  frequency  of  the  uniform  radial 
mode.   Consequently,  these  modes  will  not  respond 


^Mode  6  is  the  mode  for  which  n  =  6.   In  the  sequel, 
the  modes  will  be  referred  to  in  this  manner. 
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TABLE  II 


NATURAL  FREQUENCIES,  NATURAL  PERIODS, 
AND  BUCKLING  PRESSURES  PREDICTED  BY 
EQUATIONS  OF  ANDERSON  AND  LINDBERG  [19 ]S 


n 

n 

T 

n 

p   xlO 

0 

1.0000 

6.28 

4 

0.1382 

45.48 

11.93 

5 

0.1168 

53.82 

5.45 

6 

0.1268 

49.56 

4.47 

7 

0.1561 

40.25 

5.08 

8 

0.1951 

32.20 

5.95 

9 

0.2426 

25.90 

7.27 

10 

0.2968 

21.17 

8.81 

11 

0.3569 

17.60 

10.53 

12 

0.4231 

14.85 

12.43 

13 

0.4952 

12.69 

14.51 

14 

0.5731 

10.96 

16.75 

15 

0.6567 

9.56 

19.17 

aAll  quantities  are  nondimensional . 
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appreciably  to  the  harmonic  oscillation  of  the  axisymmetric 

2    .      . 
motion.    Discarding  the  harmonic  terms  gives 

2   --L 

0  ,  .    ^o  o     o 
w  (t)  =  j   e 

1  +  T 

o 


Since  t  >>1  for  the  loadings  of  interest  to  the  present 
study,  the  above  equation  reduces  to 

T 

0,  ,        To 
w  (t)  =  p  e 

^o 


which  represents  the  quasi-static  response  of  the  axi- 
symmetric mode  to  the  applied  load.  Substituting  this 
result  into  (5.7)  yields 

T 

=  n   pQe  °         (5.11) 

Wl  Wl 


As  pointed  out  earlier,  the  solution  of  (5.11)  is 

-t/t 

on 
hyperbolic  in  nature  when  p  e      >p     Hence,  as  long  as 

-t/t 

the  pressure  p  e      exceeds  the  critical  pressure  of  any 

one  of  the  modes,  that  mode  will  grow  in  an  exponential 

manner.   Such  modes  will  be  referred  to  as  "hyperbolic" 

modes . 

By  virtue  of  their  simplified  analysis,  Anderson  and 

Lindberg  [19]  were  able  to  express  the  response  of  the  n 

mode  in  terms  of  the  amplification  ratio  w  /wi  ,  as  can  be 


T 

••n              „ 
w                2  , 

+    n     (p 

.n                rcr 

T 
O, 

•     P0G            } 

n 
w 

.  n 

2 

The  results  of  Mclvor  and  Lovell  [20]  show  this  to 

be  true . 
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seen  in  (5.6).   They  showed  that  an  amplification  of  1000 
is  tantamont  to  permanent  buckling.   The  equations  of 
the  present  analysis  can  not  be  solved  directly  for  the 
amplification  due  to  the  nonlinear  terms  in  the  equations. 
However,  once  a  solution  is  obtained,  the  amplification 
of  each  mode  can  be  computed,  based  upon  the  values  given 
to  wi  ;  and  direct  comparisons  can  be  made  with  the 
results  of  Reference  [19]. 

Mathieu  Instability.   A  mode  can  also  be  excited  to 

large  amplitude  through  the  mechanism  of  parametric 

resonance,  or  Mathieu  instability.   For  x  >>1,  the  essen- 

2  o 

tial  features  of  the  mechanism  are  present  when  (5.10)  is 

reduced  to 

0  -t/t 

w  (t)  =  p   (e       -  cos  t)  (5.12) 

Substituting  (5.12)  into  (5.7)  and  making  use  of  (5.9)  lead 

to 

—  t/t 
••n         /r>2        2  o.2  \n  /it-ion 

w      +    (ft      -n   p   e  +np      cos    t)    w  (5.13) 

n  o  o 

-,  -t/t 

2         ,        '    o  N       .  n 

n   p       (e  -    cos    t)    wi 

o 

3 
If  t>>t   and  wi   =0,   (5.13)  reduces  to  Mathieu' s  equation 

wn  +  (ft2  +  n2p   cos  t)  wn  =  0  (5.14) 

n      o 


^The  assumption  that  t>>t   is  not  valid  for  the 
present  problem.   Consequently,  the  stability  chart  of 
(5.14)  does  not  adequately  identify  the  Mathieu  modes. 
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There  are  two  possible  solutions  to  (5.14),  both  of  them 

oscillatory.   However,  one  of  them  grows  unbounded  with 

increasing  time,  while  the  other  remains  bounded.   The 

unbounded  growth  is  known  as  parametric  resonance  or 

dynamic  instability.   The  name  Mathieu  mode,  as  used  in 

the  sequel,  is  given  to  any  mode  that  grows  to  large 

amplitude  in  an  oscillatory  manner. 

The  stability  boundaries  for  the  solutions  of  (5.14) 

are  well  documented  [7,  47]  and  are  shown  by  Figure  5.1. 

2    2 
If  a  point  (ft  ,  n  p  )  falls  in  the  unstable  region,  the 

solution  for  that  mode  will  grow  without  bound.   The 
stability  chart  shows  that  when  the  frequency  of  a 
flexural  mode  is  approximately  one-half  of  the  axisymmetric 
natural  frequency  (ft-  =  1.0),  the  response  of  that  flexural 
mode  will  be  unstable,  even  for  a  very  small  peak  pressure. 
The  points  corresponding  to  modes  12,  13,  14,  and  15  for 
the  shell  considered  here  are  shown  in  Figure  5.1. 

Anderson  and  Lindberg  [19]  were  interested  in  pre- 
dicting the  critical  buckling  loads.   From  experimental 
results,  they  knew  that  dynamic  buckling  occurs  as  the 
result  of  very  rapid  growth  of  the  buckling  mode  accom- 
panied by  very  little  oscillatory  motion.   However,  the 
Mathieu  instability  occurs  in  an  oscillatory  manner  re- 
quiring several  cycles  before  the  displacements  become 
large.   Furthermore,  when  the  amplitude  of  a  Mathieu  mode 
becomes  large  enough  to  cause  plastic  flow,  energy  will 
be  dissipated,  the  axisymmetric  oscillations  will  decay, 
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FIGURE  5. 1 

MATHIEU  STABLITY  CHART  FOR  HARMONIC 
UNIFORM  PRESSURE 
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and  subsequent  growth  of  the  Mathieu  mode  would  most 

4 
likely  not  occur.    Consequently,  they  were  able  to 

disregard  the  response  of  the  Mathieu  modes. 

Mclvor  and  Lovell's  Contribution 

Under  impulsive  loads,  the  initial  growth  of  the 
flexural  modes  is  governed  by  the  Mathieu  equation.   The 
impulsive  response  of  the  axisymmetric  mode  can  be  ob- 
tained from  (5.10)  by  taking  the  limit  as  t   goes  to  zero 
while  the  product  p  x   remains  constant.   The  result  is 

w°  (T)  =  P0TQ  sin  x  (5.15) 

Substituting  (5.15)  into  (5.7)  and  making  use  of  (5.9) 
lead  to 

••n  ,  ,n2  2  „.  x  n     .n  2  ,_  ,,s 

w   +  {u      -  n  p  x   sin  x)w   =  wi  n  p  x   sin  x   (5.16) 
n     ^o  o  *o  o 

If  wi   =0,  (5.16)  reduces  to  the  Mathieu  equation 


wn  +  (ft2  -  n2p  x   sin  x)wn  =  0  (5.17) 


Mclvor  and  Lovell's  investigation  [20]  was  concerned 
with  the  undamped  response  to  purely  impulsive  loads.   With 
an  equation  similar  to  (5.17),  they  identified  the  Mathieu 
modes  from  a  stability  chart.   The  total  deflection  was 
assumed   to  consist  of  the  Mathieu  modes  plus  the  axi- 
symmetric mode.   An  energy  balance  revealed  that  the 


The  axisymmetric  oscillations  are  due  to  the  initial 
conditions.   Thus,  there  is  a  limited  amount  of  energy  in 
them.   Once  this  energy  is  dissipated,  they  will  no  longer 
exist;  and  the  mechanism  for  dynamic  instability  is  removed 
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response  of  the  neglected  modes  was,  in  fact,  negligible. 
All  of  the  nonlinear  coupling  between  the  modes  was  re- 
tained.  During  the  response,  the  modes  underwent  beat-like 
oscillations  accompanied  by  a  cyclic  exchange  of  energy 
between  the  modes.   The  maximum  amplitude  of  a  mode  did 
not  always  occur  during  the  first  beat,  sometimes  occurring 
during  the  second  beat.   The  stresses  developed  during  the 
response  to  an  impulse  were  several  times  higher  than  the 
stresses  of  an  axisymmetric  response  to  the  same  impulse. 
Furthermore,  the  largest  stress  usually  did  not  occur 
during  the  first  significant  exchange  of  energy,  but 
occurred  later  when  several  modes  were  nearly  in  phase. 
The  maximum  amplitudes  attained  by  the  flexural  modes  were 
nearly  independent  of  the  magnitudes  of  the  perturbation 
in  each  mode;  however,  the  phasing  between  the  modes  was 
affected  by  the  magnitudes  of  the  perturbation. 

III.   INVESTIGATION  OF  THE  PROBLEM  BY  THE 
METHODS  DEVELOPED  IN  THIS  STUDY 

Preliminaries 

Minimum  Number  of  Axial  Grid  Points.   In  order  to 
minimize  the  computing  time,  the  number  of  axial  stations 
was  kept  as  small  as  reasonable  accuracy  permitted.   The 
criteria  used  for  assessing  the  adequacy  of  a  given  number 
of  stations  were:  (1)  the  accuracy  of  the  circumferential 
in-plane  forces,  (2)  the  accuracy  of  the  amplification 
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factors,  and  (3)  the  accuracy  of  the  representation  of  the 
axial  variation  of  the  radial  displacements. 

Since  the  load  and  the  initial  imperfection  were 
assumed  to  be  symmetric  with  respect  to  the  midspan  of  the 
shell,  it  would  not  be  unreasonable  to  expect  the  response 
to  also  be  symmetric  with  respect  to  the  midspan.   If  true, 
this  would  effect  a  considerable  savings  in  the  computing 
time.   Two  fully  nonlinear  solutions  were  computed  to 
determine  whether  or  not  the  response  was  indeed  symmetric. 
The  first  run  employed  the  simple  support  conditons  at 
both  ends  of  the  shell  and  fourteen  axial  stations.   The 
second  run  used  the  midspan  symmetry  conditons  and  seven 
axial  stations.   The  displacements  of  the  two  solutions 
were  identical  to  at  least  four  significant  figures. 
Therefore,  the  response  for  this  problem  was  assumed  to  be 
symmetric  with  respect  to  the  midspan  of  the  shell. 

A  series  of  computer  runs  was  made  using  a  different 
number  of  axial  stations  over  the  half-length  of  the  shell 
for  each  run.   There  were  no  significant  differences  in 
the  results  when  nine,  eleven,  or  thirteen  axial  stations 
were  used  along  one-half  of  the  shell.   However,  when  the 
number  of  stations  was  reduced  to  seven,  the  error  in  the 
circumferential  in-plane  force  was  approximately  7  per  cent 
near  the  end  of  shell;  and  smaller  errors  were  noted  in 
the  radial  displacements  and  the  amplification  factors. 
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Consequently,  nine  stations  over  the  half-span  of  the  shell 
were  used  for  the  solutions  presented  in  the  following 
sections  of  this  chapter. 

Selection  of  the  Participating  Modes.   Since  the 

computer  program  is  dimensioned  to  accomodate  a  maximum 

5 
of  fifteen  Fourier  modes,   the  series  expansions  for  the 

displacements  had  to  be  truncated.   However,  the  truncated 

series  was  not  restricted  to  the  first  fifteen  modes, 

modes  0  through  14.   Instead,  the  participating  modes 

were  selected  to  include  only  those  expected  to  respond 

to  the  particular  loading  condition.   The  procedure  used 

for  selecting  the  responding  modes  can  best  be  explained 

by  an  example. 

Consider  the  loading 

P(t)  =  70  e"t/,0°3  (5.18) 

-4 
Comparing  the  70  psi  (p   =  6.37  x  10   )  peak  pressure  with 

the  buckling  pressures  listed  in  Table  II  reveals  that  modes 
5  through  8  should  exhibit  hyperbolic  behavior  during  the 
early  phases  of  the  response.   Hence,  these  modes  were 
included  in  the  series . 

Identifying  the  Mathieu  modes  is  a  little  more  diffi- 
cult.  As  a  result  of  the  assumptions  made  to  arrive  at 
(5.14),  the  stability  chart  of  Figure  5.1  is  not  sufficient 


5 
The  number  of  modes  is  limited  primarily  by  com- 


puting time 
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to  identify  the  Mathieu  modes.   Consequently,  a  procedure 
was  used  that  investigated  the  stability  of  each  of  the 
possible  Mathieu  modes.   For  example,  Table  II  shows  that 
modes  12  through  15  are  possible  Mathieu  modes  since  they 
have  frequencies  lying  in  the  neighborhood  of  one-half 
the  frequency  of  the  axisymmetric  mode.   To  determine 
which  of  these  modes  were  Mathieu  modes,  the  series  was 
truncated  so  that  only  modes  0  and  12  through  15  were 
retained.   As  a  result  of  this  choice  of  Fourier  indices, 
each  of  the  flexural  modes  couples  only  with  mode  zero. 
Finally,  the  coupling  terms  in  the  equations  of  mode  0 
were  equated  to  zero  by  a  simple  alteration  of  the  computer 
program.   This  eliminated  the  mechanism  for  the  extraction 
of  energy  from  the  mode  zero  oscillation.   The  resulting 
equations  of  motion  of  the  flexural  modes  resembled 
Mathieu 's  equation,  and  any  mode  that  grew  to  a  large 
amplitude  under  these  circumstances  was  assumed  to  be  a 
Mathieu  mode.   A  computer  run  was  carried  out  based  on 
this  procedure  using  the  pressure  given  by  (5.18).   The 


The  modes  which  couple  are  determined  by  forming  all 
combinations  of  sums  and  differences  of  the  Fourier  numbers 
in  the  series.   If  any  sum  or  difference  equals  any  of  the 
Fourier  numbers  in  the  series,  then  the  two  Fourier  numbers 
contributing  to  that  sum  or  difference  are  said  to  couple 
with  the  third  mode.   For  example,  modes  3  and  5  couple  to 
influence  modes  2  and  8  since  3+5=8  and  5-3=2.   Mode 
13  couples  with  itself  to  influence  modes  0  and  26  since 
13  -  13  =  0  and  13  +  13  =  26.   For  a  more  detailed  explana- 
tion see  Appendix  B,  Reference  [30]. 
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results  are  given  in  Table  III.  The  maximum  amplification 
listed  in  the  table  is  defined  as  the  maximum  value  of  the 
ratio 


n  i 

|  w .  | 

—  (5.19) 

— n  i 

w. 

1  * 


where  station  i  is  the  whole  station  nearest  to  the  midspan 
of  the  shell.   From  Table  III,  mode  14  can  be  readily 
identified  as  a  Mathieu  mode.   However,  as  shown  in  Figure 
5.1,  the  Mathieu  stability  chart  of  (5.14)  indicates  that 
mode  13  should  be  the  unstable  mode. 

A  possible  explanation  of  this  discrepancy  can  be 
surmised  by  considering  (5.13).   When  wi   is  equal  to  zero, 
(5.13)  reduces  to 

..n    /r.2     2        o     2         .  n 

w   +  (ft   -  n  p  e       +np   cos  t)w   =  0      (5.20) 
n      ^o  ^o 

TABLE  III 
RESPONSE  OF  INDIVIDUAL  MODES  WITHOUT 
EXTRACTION  OF  ENERGY  FROM  AXISYMMETRIC  MODE 

p   =  6.37  x  10~4,  T   =  208 
*o  o 

Maximum  Amplification 
Mode     Attained  During  Interval 
x  =  0  to  140 

12  1.8 

13  1.8 

14  1990.0 

15  3.4 
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This  equation  is  similar  to  the  Mathieu  equation,  except 

2 

that  ft   has  been  replaced  by 

n2     2    "t/to  (5.21) 

ft   -  n  p  e 
n      ^o 

The  quantity  expressed  by  (5.21)  can  be  thought  of  as  a 
time  varying  frequency  of  the  n  '  mode.   If  this  quantity 

is  plotted  along  the  abscissa  of  the  Mathieu  stability 

2 

chart  instead  of  ft  ,  as  in  Figure  5.2,  the  points  associated 

with  the  flexural  modes  12,  13,  14  and  15  move  horizontally 
across  the  chart  as  r  increases. 

Examination  of  Figure  5.2  reveals  that  at  t  =  0 , 
modes  12  and  13  are  in  the  stable  regions,  and  modes  14  and 
15  are  unstable.   In  computing  the  amplifications  listed 
in  Table  III,  the  integration  was  terminated  at  t  =  140. 
During  that  time  interval,  modes  12  and  13  grew  very  little; 
mode  14  was  amplified  nearly  two  thousand  times;  and  mode 
15  was  amplified  3.4  times.   These  results  can  be  correlated 
with  the  trajectories  of  the  corresponding  points  in  Figure 
5.2.   During  the  interval  of  integration,  modes  12  and  13 
remained  in  a  stable  region,  but  mode  14  remained  in  the 
unstable  region.   Mode  15,  which  was  initially  in  an  un- 
stable region,  quickly  moved  into  a  stable  region.   Since 
growth  to  a  large  amplitude  requires  several  oscillations 
of  the  axisymmetric  mode,  mode  15  moved  into  a  stable  region 
before  having  time  to  grow  appreciably.   Had  the  integration 
continued,  mode  13  would  have  entered  the  unstable  region; 
but  mode  12  would  have  remained  stable. 
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FIGURE  5.2 
MATHIEU  STABILITY  CHART  FOR  EXPONENTIALLY 
DECAYING  UNIFORM  PRESSURE 

P  =70  pti(p  =6.3  7xl0~4),  t=  3.0msec  (t  =  208) 

0  0  0  0 


128 


As  a  consequence  of  this  correlation  the  stability 
chart  of  Figure  5.2,  which  is  based  on  the  simplified 
equations  of  Reference  [19],  was  considered  to  be  adequate 
for  identifying  the  Mathieu  modes  for  the  problem  under 
consideration.   To  identify  the  Mathieu  modes  for  a  prob- 
lem with  different  loading  conditions,  it  might  be  neces- 
sary to  make  several  preliminary  computer  runs  similar  to 
the  one  described  above. 

When  computing  a  fully  coupled  nonlinear  solution 
for  the  loading  given  by  (5.18),  modes  0,  1,  4,  5,  6,  7, 
8,  9,  10,  13,  14,  and  15  were  retained  in  the  Fourier 
expansions.   Note  that  modes  4,  9,  and  10,  modes  in  the 
neighborhood  of  the  hyperbolic  modes,  were  retained.   The 
Mathieu  modes  that  were  retained  included  all  modes  whose 
trajectories  in  Figure  5.2  passed  through  the  unstable 
region.   Finally,  mode  1  was  also  included  because  it  be- 
comes involved  in  the  coupling  of  all  the  modes  whose 
wave  numbers  differ  by  unity.    The  adequacy  of  this  choice 
of  modes  was  established  by  computing  a  solution  in  which 
modes  0  through  14  were  retained.   The  extra  modes,  2,  3, 
11,  and  12,  did  not  respond  appreciably;  and  including  them 
did  not  alter  the  response  of  the  other  modes. 


7 
In  the  Donnell  theory  mode  1  is  inaccurately 

described  for  many  classes  of  shells.   However,  if  the  shell 

is  thin  and  L/as2 ,  as  is  the  case  here,  the  representation 

of  mode  1  is  reasonably  accurate  [48]. 
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Results 

Computer  solutions  were  generated  for  four  specific 

o 

cases,  hereafter  referred  to  as  A,  B,  C,  and  D.    In  cases 

A  and  B,  the  damping  coefficient  c  was  taken  equal  to  zero, 

-5 
while  in  cases  C  and  D,  c  =  3.416  x  10   .   Cases  A  and  C 

were  each  composed  of  three  runs  with  t   =3.0  msec 

(t   =  208)  and  P   equal  to  60  psi  (p   =  5.46  x  10~4),  70 
o  o  c  ^o 

psi  (p   =  6.37  x  10~4) ,  and  77  psi  (p   =  7.01  x  10~4). 

Cases  B  and  D  were  also  composed  of  three  runs  each,  but 

at  three  different  peak  pressures  with  t   =1.0  msec 

^     ^  o 

(t   =  69.3).   The  results  for  these  cases  will  be  compared 
with  the  results  obtained  by  Anderson  and  Lindberg  [19] 
and  Mclvor  and  Lovell  [20]. 

The  significance  of  the  amplitudes  of  the  initial 
imperfection  modes  was  also  investigated  by  computing  a 
solution  in  which  the  amplitudes  wi   were  taken  from  a 
table  of  random  numbers.   Finally,  some  results  are  pre- 
sented to  show  the  importance  of  the  nonlinear  coupling 
between  the  modes 


Case  A.   The  results  for  this  series  of  runs  are 
summarized  in  Table  IV.   The  first  column  shows  the  Fourier 
indices  of  the  modes  included  in  the  solution.   The  second 
column  shows  the  maximum  amplification  of  each  flexural 
mode,  as  given  by  (5.19).   The  third  column  contains  the 


A  typical  solution  required  approximately  45  minutes 
on  the  IBM  360/67. 


130 


TABLE  IV 


CASK  A,  AMPLIFICATIONS  AND  STRESSES 

t   =  3.0  msec  (x   =  208), ZERO  DAMPING 
o  o 


P   -  60 
o 

psi 

:  po  = 

=  5.. 

46  x 

10 

•4 

Maximum 

circumferential 

stress 

=  21, 

r600 

psi 

Maximum 

circumferential 

stress 

occurred 

at  t 

= 

116 

Integration 

terminated  at  t  = 

118 

mode 

maximum 

amp . 

T 

for 

maximum  amp. 

0 

1  .72  110.4 

4  1.53  27.83 

5  11.22  52.67 

6  41.02  72.05 

7  24.15  52.79 

8  9.18  33.77 

9  4.60  20.13 

13  1.94  117.6 

14  97.39  116.4 

15  1.56  116.6 


P   =  70  psi  :  p   =  6.37  x  10  4 

o       r  ^o 


Maximum  circumferential  stress  =  32,200 

psi 

Maximum  circumferential  stress  occurred 

at  t  =  165 

Integration  terminated  at  x  =  198 

mode       maximum  amp .       t  for 

maximum  amp . 

0 

1  7.93  190.6 

4  2.22  88.35 

5  26.70  71.71 

6  222.3  97.35 

7  101.7  72.23 

8  22.02  40.37 

9  7.84                  26.41 
10  4.15                  19.43 

13  3.80  193.8 

14  229.2  164.6 

15  4.06  162.8 
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TABLE  IV  Continues 


P   =  77 
o 

psi 

:  p 

=  7.i 

31  x 

10~ 

■4 

Maximum 

circumferential 

stress 

=  64, 

300 

psi 

Maximum 

circumferential 

stress 

occurred 

at  t 

= 

125 

Integrat 

:ion 

terminated  < 

at  t  = 

166 

mc 

>de 

maximum 

amp . 

T 

for 

maximum  amp . 

0 

1  126.8  116.7 

4  10.90  159.7 

5  134.7  135.0 

6  1411.  127.4 

7  423.8  90.96 

8  97.20  118.3 

9  25.80  118.9 
10  7.79  145.3 

13  116.6  145.4 

14  239.5  144.1 

15  41.17  157.7 


nondimensional  time  at  which  the  maximum  amplification 
occurred.   The  magnitude  of  the  maximum  circumferential 
stress  at  6  =  0  ,  as  given  by  (4.31) ,  and  the  nondimen- 
sional time  of  its  occurrence  are  also  listed.   Finally, 
the  nondimensional  time  at  which  the  run  was  terminated 
is  indicated. 

The  time  histories  of  the  radial  displacements  at 
the  whole  station  nearest  to  the  midspan  are  shown  in 
Figures  5.3-5.5  for  the  four  or  five  modes  that  attained 
the  largest  amplifications.   In  addition,  the  Figures  are 


9 
The  nondimensional  period  of  the  axisymmetric  mode 

is  approximately  equal  to  2 tt  . 
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annotated  with  values  of  the  circumferential  stress  to 
show  the  times  at  which  a  relative  maximum  occurred.   Note 
the  different  scales  for  w   in  the  three  Figures . 

From  Table  IV,  it  can  be  seen  that  when  the  peak 

-4 
pressure  was  60  psi  (p   =  5.46  x  10   ) ,  hyperbolic  modes 

5,  6,  and  7  grew  to  large  amplitudes,  and  that  the  growth 
occurred  during  the  early  phases  of  the  response.   Exam- 
ination of  Table  II  shows  that  these  are  precisely  the 
modes  having  static  buckling  pressures  less  than  the  peak 
pressure  of  the  load  pulse.   Furthermore,  mode  6,  the 
static  buckling  mode,  attained  a  larger  amplification  than 
the  other  hyperbolic  modes.   If  comparisons  are  made  for 
the  other  two  runs,  similar  observations  can  be  made. 

In  the  third  run,  mode  1  grew  to  an  amplification  of 
126.8  just  prior  to  the  peaking  out  of  the  hyperbolic  modes. 
This  growth  of  mode  1  was  very  likely  due  to  coupling 
between  the  hyperbolic  modes  that  produced  large  coupling 

terms  in  the  mode  1  equations.   The  reason  for  including 
mode  1  in  the  truncated  series  is  now  clear. 

In  all  three  runs,  mode  14  attained  the  largest  am- 
plification of  the  Mathieu  modes.   Inspection  of  Figures 
5.3-5.5  reveals  that  mode  14  underwent  beat-like  oscilla- 
tions requiring  several  cycles  to  attain  a  maximum  ampli- 
tude.  However,  it  was  not  clear  in  any  run  that  the  largest 

amplitude  of  the  Mathieu  modes  had  occurred  prior  to  the 
termination  of  the  computer  run.   In  fact,  the  information 
shown  in  Figure  5.2  indicates  that  mode  13  would  have  grown 
large  for  x>i     . 
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Comparison  of  the  time  histories  for  the  three  runs 
reveals  that  at  the  lower  peak  pressures,  the  Mathieu 
mode  dominated  the  response.   As  shown  in  Figure  5.3,  the 
amplification  attained  by  mode  14  was  more  than  twice  as 
large  as  the  amplification  of  any  of  the  hyperbolic  modes. 
Furthermore,  the  largest  stresses  occurred  as  a  result  of 
the  growth  of  the  Mathieu  mode.   As  the  peak  pressure  and 
the  total  impulse  P  t   were  increased,  the  amplifications 
of  both  the  Mathieu  mode  and  hyperbolic  modes  were  of 
about  equal  magnitude,  as  shown  in  Figure  5.4.   Figure 
5.5  shows  that  at  the  highest  pressure  and  impulse,  the 
hyperbolic  modes  dominated,  and  the  maximum  stresses 
occurred  as  a  result  of  their  growth. 

In  Figure  5.6,  the  maximum  amplification  of  the  hyper- 
bolic mode  6  is  plotted  against  the  peak  pressure  and, 
equivalently ,  the  total  impulse.   Assuming  an  amplification 
of  1000  to  be  the  dynamic  buckling  limit,    Figure  5.6 
shows  that  the  buckling  pressure  is  76  psi  when  the  total 
impulse  is  228  psi-msec.   Note  the  logarithmic  scale  for 
the  amplification  factor. 

Case  B.   For  the  runs  of  case  B,  the  time  constant 

of  the  decaying  load  was  reduced  to  1.0  msec  (t   =  69.3); 

and  the  pressures  for  the  three  runs  were  100  psi 

(p   =  9.10  x  10~4;,  110  psi  (p   =  10.01  x  10~4),  and  120 

-4 
psi  (p   =  10.92  x  10   ).   The  hyperbolic  modes  were  5,  6, 


Anderson  and  Lindberg  [19]  also  used  this  criterion 
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7 ,  8,  and  9,  all  of  which  have  buckling  pressures  less 
than  the  peak  pressures  of  the  loads.   In  addition,  modes 
4  and  10  were  included.   The   Mathieu   modes  were  iden- 
tified from  a  stability  chart  similar  to  Figure  5.2  and 
included  modes  13,  14,  and  15  at  the  lowest  peak  pressure 
plus  mode  12  at  the  two  higher  peak  pressures.   As  with 
case  A,  the  damping  coefficient  was  set  to  zero.   Two  of 
the  runs  were  computed  out  to  T    =  20  0;  and  the  third,  to 
t  =  150.   The  results  are  summarized  in  Table  V. 

As  in  Case  A,  the  character  of  the  response  changed 

considerably  as  the  peak  pressure  and  the  total  impulse 

-4 
were  increased.   When  P   =  100  psi  (p   =  9.10  x  10   ) , 

o       ^    *o 

the   maximum  amplifications  of  the  Mathieu  modes  were 
roughly  twice  as  large  as  the  amplifications  of  the 
hyperbolic  modes.   For  the  second  solution,  the  peak  pres- 
sure was  equal  to  110  psi;  and  the  maximum  amplifications 
of  the  hyperbolic  and  Mathieu  modes  were  approximately 
equal.   The  loading  of  the  final  run  was  near  the  critical 
dynamic  buckling  limit  and  boosted  the  hyperbolic  mode  6 
to  an  amplification  of  832.6.   The  largest  amplification 
of  the  Mathieu  modes  was  340.6.   The  maximum  stress  occurred 
well  after  the  hyperbolic  modes  had  passed  their  peak 
amplifications  in  the  first  and  second  runs.   On  the  other 
hand,  the  largest  stress  occurred  in  conjunction  with  the 
large  response  of  the  hyperbolic  modes  in  the  last  run. 
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TABLE  V 

CASE  B,  AMPLIFICATIONS  AND  STRESSES 

t   =1.0  msec  (T   =69.3),  ZERO  DAMPING 
o  o 


P   =  100  psi 
o        ^ 


o 


9.10  x  10 


-4 


Maximum  circumferential  stress  =  55,000  psi 


Maximum  circumferential  stress  occurred  at  t 


199 


Integration  terminated  at  t  =  200 


mode 


maximum  amp 


t  for  maximum  amp 


0 

1 

4 

5 

6 

7 

8 

9 

10 

11 

13 

14 

15 


96 
5 

32 
134 
119 

51 

20 
9 

13 
331 
245 

73 


15 

34 

85 

2 

7 

7  7 

38 

64 

85 

0 

1 

04 


197 

185 

58 

65 

52 

40 

156 

19 

182 

196 

77 

193 


7 

6 

55 

44 

72 

24 

6 

79 

9 

2 

09 


P   =  110  psi                    :  p   =10, 
o        ^                       ^o 

.01  x  10  4 

Maximum  circumferential  stress  =  65,500 

psi 

Maximum  circumferential  stress  occurred 

at  t  =  117 

Integration  terminated  at  t  =  200 

mode       maximum  amp.       t  for 

maximum  amp. 

0 

1 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 


77.55 
14.82 
57.19 
291.5 
256.7 
106.1 
54.66 
17.88 
51.86 
226.5 
273.2 
258.7 
138.6 


125 
179 


63.56 
71.22 


59 
46 
172 
131 
195 
184 
184 
112 
169 


28 

52 

8 

4 

6 

2 

1 

5 

6 
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TABLE  V   Continued 


r   =  120  psi                   :  p   =10 

o        v                                                       ^o 

.92  x  10"4 

Maximum  circumferential  stress  =  76,200 

psi 

Maximum  circumferential  stress  occurred 

at  t  =  57 

Integration  terminated  at  t  =  150 

mode       maximum  amp.       t  for 

maximum  amp. 

0 

1  148.6  78.31 

3  36.77  100.3 

4  34.58  136.3 

5  160.5  81.54 

6  832.6  82.50 

7  656.4  66.93 

8  218.2  52.06 

9  70.54  77.54 

10  45.48  128.4 

11  30.94  135.7 

12  152.8  114.2 

13  274.4  141.0 

14  340.6  113.6 

15  250.3  98.83 


The  behavior  of  the  Mathieu  modes  reflected  the 
shorter  time  constant  and  the  higher  peak  pressures.   In 
all  three  runs,  all  of  the  included  Mathieu  modes  grew  to 
significant   amplitudes.   Again  this  can  be  explained  by 
the  stability  chart  of  Figure  5.2.   As  n  p   is  increased, 
the  trajectories  of  the  Mathieu  modes  are  shifted  upward 
and  extended  in  length  to  the  left;  the  reduced  time  con- 
stant causes  the  points  corresponding  to  the  modes  to 
move  across  the  chart  at  a  faster  rate.   The  net  result 
was  that  modes  12,  13,  and  15,  as  well  as  14,  were  in  the 
unstable  region  sufficiently  long  to  respond. 
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Note  that  in  all  three  runs,  the  stress  exceeded  42 
ksi,  the  yield  point  stress  of  the  shell  material,  while 
the  amplifications  of  the  hyperbolic  modes  never  exceeded 
1000,  the  dynamic  buckling  limit  used  by  Anderson  and  Lind- 
berg  [19]. 

The  maximum  amplification  of  the  hyperbolic  modes  as 
a  function  of  the  peak  pressure  and  the  total  impulse  is 
shown  in  Figure  5.7.   An  amplification  of  1000  corresponds 
to  a  pressure  of  approximately  122  psi. 

Comparison  with  Anderson  and  Lindberg's  Findings. 

Anderson  and  Lindberg  [19]  identified  the  peak  pressure 

P   and  the  total  impulse  P  t   as  the  significant  parameters 
o  o  o 

in  the  dynamic  buckling  problem.   Their  dynamic  buckling 
curve,  along  which  the  maximum  amplification  of  the  hyper- 
bolic modes  was  1000,  is  shown  in  Figure  5.8.   The  sharp 
break  in  the  curve  divides  the  regime  of  the  elastic  model, 
on  the  lower  right,  from  the  regime  of  the  tangent  modulus 
model,  at  the  upper  left.   As  the  impulse  of  the  load  goes 
to  infinity,  corresponding  to  a  step  load,  the  dynamic 
buckling  curve  approaches  the  static  buckling  pressure 
line.  The  dynamic  curve  shows  that  the  shell  can  endure 
a  pulse  having  a  peak  pressure  larger  than  the-  static 
buckling  pressure,  with  the  critical  peak  pressure  depend- 
ing on  the  total  impulse  of  the  load.   While  Anderson  and 
Lindberg's  model  gives  the  dynamic  buckling  curve,  it  does 
not  predict  the  response  in  the  region  lying  below  the 
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dynamic  buckling  curve  since  the  Mathieu  modes  were  not 
included.   However,  the  algorithm  of  this  dissertation  is 
capable  of  computing  the  displacements  and  the  stresses 
in  the  subcritical  region,  notably,  in  the  region  between 
the  static  buckling  line  and  the  dynamic  buckling  curve. 

The  two  straight  lines  running  diagonally  across 
the  Figure  are  lines  along  which  t   is  equal  to  3.0  msec, 
case  A,  and  1.0  msec,  case  B.   The  two  points  labeled 
CASE  A  and  CASE  B  are  the  buckling  pressures  and  impulses 
which  were  obtained  from  Figures  5.6  and  5.7.   They  are 
based  on  an  amplification  of  1000.   Note  that  the  points 
fall  very  close  to  Anderson  and  Lindberg's  buckling  cruve . 
Therefore,  Anderson  and  Lindberg's  treatment  of  the  boundary 
conditions  and  their  neglecting  of  certain  nonlinear 
coupling  terms  do  not  lead  to  any  signigicant  errors  in 
predicting  the  amplification  of  1000. 

Comparison  with  Mclvor  and  Lovell's  Findings.   Under 
an  impulsive  load,  as  considered  by  Mclvor  and  Lovell  [20], 
the  axisymmetric  mode  oscillates  about  its  unloaded  static 
equilibrium  state,  and  the  quasi-static  response  of  the 
axisymmetric  mode  is  zero.   Therefore,  the  hyperbolic 
modes,  which  respond  primarily  to  the  quasi-static  motion 
of  the  axisymmetric,  are  not  excited.   In  cases  A  and  B, 
the  amplitudes  of  the  hyperbolic  modes  were  small  at  the 
lower  pressures  and  impulses,  because  the  quasi-static 
response  of  the  asisymmetric  mode  was  very  small;  yet  the 
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high  frequency  oscillation  of  the  axi symmetric  mode  was 
sufficient  to  excite  the  Mathieu  modes.   The  observed 

behavior  of  the  Mathieu  modes  agrees  qualitively  with 

11 
Mclvor  and  Lovell's  findings.     For  instance,  in  Figure 

5.4,  the  Mathieu  mode  responded  in  a  beat-like  manner 
requiring  several  cycles  to  attain  the  first  peak.   Fur- 
thermore, the  peak  amplitude  of  the  Mathieu  mode  was 
larger  in  the  second  beat  than  in  the  first,  leading  to  a 
large  stress. 

Maximum  Stresses.   Although  the  theory  of  Reference 
[19]  appears  to  be  very  adequate  for  computing  buckling 
thresholds,  the  results  for  cases  A  and  B,  contained  in 
Tables  IV  and  V,  indicate  that  the  stresses  may  exceed 
the  yield  point  stress  at  peak  pressures  and  total  impulses 
well  below  the  buckling  threshold.   When  the  load  was  sub- 
critical,  maximum  stress  occurred  well  after  the  lower 
frequency,  buckling  modes  had  attained  their  maximum 
amplitudes.   However,  at  loads  near  critical  the  maximum 
stress  occurred  at  about  the  time  the  dominant  hyperbolic 
mode  was  at  its  largest  amplitude. 

Since  a  structure  designed  for  repeated  use,  such 
as  the  NWL  shock  tube ,  would  experience  fatigue  weakening 
if  the  stresses  exceed  the  yield  point  stress,  one  might 
ask  whether  a  normal  amount  of  viscous  damping  would 


A  quantitative  comparison  can  not  be  made  due  to 
the  differences  in  the  analyses. 
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attenuate  the  Mathieu  response,  thus  reducing  the  stresses. 
This  aspect  is  considered  in  cases  C  and  D. 

Case  C.   The  loading  of  case  C  was  the  same  as  the 
loading  of  case  A  except  that  viscous  damping  was  included. 

The  nondimensional  coefficient  c  was  taken  equal  to 

-5 
3.416  x  10   .   The  maximum  amplification  of  each  mode  and 

the  time  of  its  occurrence,  the  maximum  circumferential 
stress  and  the  time  of  its  occurrence,  and  the  time  at  which 
the  integration  was  terminated  are  contained  in  Table  VI. 
In  addition,  Figures  5.9-5.11  show  the  time  histories  of 
the  active  modes,  annotated  to  indicate  the  times  at  which 
large  stresses  occurred. 

The  computer  runs  of  case  C  were  carried  out  for 
longer  periods  of  time  than  were  the  runs  of  case  A  in 
order  to  include  any  significant  growth  of  mode  13  that 
might  lead  to  high  stresses.   The  amplifications  of 
Table  VI  indicate  that  mode  13  had  moved  into  the  unstable 
region  of  the  Mathieu  chart  of  Figure  5.2  prior  to  the  end 
of  each  run.   On  the  first  two  runs,  mode  13  became  un- 
stable only  after  the  damping  effects  had  extracted  much 
of  the  energy  from  the  motion  of  the  shell.   Consequently, 
mode  13  did  not  grow  sufficiently  to  generate  a  higher 
stress  than  had  previously  occurred.   Due  to  the  higher 
peak  pressure  of  the  third  run,  mode  13  became  unstable 
much  sooner.   However,  the  mode  13  response  was  not  respon- 
sible for  the  maximum  stress  since  the  motion  of  the 
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TABLE  VI 


CASE  C,  AMPLIFICATIONS  AND  STRESSES 

t   =  3.0  msec  (r   =  208),  DAMPING  INCLUDED 
o  o 


-4 
P   =  60  psi  :  p   =  5.46  x  10 

Maximum  circumferential  stress  =  16,250 

Maximum  circumferential  stress  occurred  at  t  =  136 

Integration  terminated  at  t  =  300 

mode      maximum  amp.       t  for  maximum  amp 

0 

1  .94  287.5 

4  1.49  27.65 

5  10.61  52.78 

6  38.24  72.01 

7  22.98  52.78 

8  8.93  33.86 

9  4.52  20.21 

13  22.86  296.0 

14  71.96  135.7 

15  1.45  84.41 


P   =  70 
o 

psi 

:  p 
^o 

=  6.3" 

'  X 

ID"4 

Maximum 

circumferential 

stress 

=  26 

,400 

psi 

Maximum 

circumferential 

stress 

occurred 

at 

T  = 

110 

Integration 

terminated  < 

at  t  = 

257 

mode 

maximum 

amp. 

T 

for 

maximum 

amp. 

0 

1 

4.53 

4 

1.99 

5 

24.79 

6 

202.0 

7 

95.82 

8 

21.26 

9 

7.68 

10 

4.09 

13 

4.98 

14 

147.5 

15 

3.09 

85.72 
33.57 
71.78 
97.66 
72.14 
46.44 
26.54 
19.41 
252.1 
177.3 
36.32 
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TABLE  VI   Continued 


P   =  77  psi  :  p   =  7.01  x  10 

o       ^  ^o 

Maximum  circumferential  stress  =  53,800  psi 

Maximum  circumferential  stress  occurred  at  r  =  110 


Integration  terminated  at  t  =  198 


0 

1 

91.06 

4 

5.33 

5 

87.00 

6 

1154. 

7 

396.2 

8 

69.72 

9 

17.71 

10 

6.54 

13 

155.8 

14 

174.1 

15 

9.92 

mode      maximum  amp.       t  for  maximum  amp 


103.9 
196.5 
119.6 
127.3 

91.54 
118.1 
119.0 
168.1 
197.0 
103.5 
183.0 


Mathieu  modes  was  greatly  overshadowed  by  the  early  growth 
of  the  hyperbolic  modes. 

The  effectiveness  of  the  damping  in  attenuating  the 
Mathieu  modes  can  be  ascertained  by  comparing  the  mode  14 
amplifications  of  the  damped  and  undamped  responses. 
Comparison  of  the  amplifications  of  mode  14  contained  in 
Table  IV  with  those  contained  in  Table  VI  reveals  that  the 
maximum  amplification  of  mode  14  was  reduced  by  2  5  to  3  0 
per  cent  by  the  inclusion  of  damping. 

Assuming  that  the  damping  would  have  extracted 
sufficient  energy  to  preclude  the  development  of  stresses 
higher  than  those  computed  in  Table  VI,  a  critical  pressure 
and  impulse  based  on  the  yield  point  stress  can  be 
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determined.   In  Figure  5.12  the  maximum  circumferential 
stress  is  shown  as  a  function  of  the  peak  pressure.   The 
critical  pressure  and  impulse  corresponding  to  a  maximum 
stress  of  42  ksi  are  75  psi  and  225  psi-msec,  respectively. 
Note  that  the  yield  stress  occurs  at  a  pressure  where  the 
hyperbolic  modes  have  much  larger  amplifications  than  the 

Mathieu  modes.   However,  the  bending  stress,  given  by  the 

2 
second  term  of  (4.29),  depends  on  n  ;  consequently,  the 

Mathieu  modes  with  their  higher  mode  numbers  can  make 

significant  contributions  to  the  stress  even  though  their 

amplitudes  may  be  small. 

Case  D.   In  case  D,  the  time  constant  t  was  taken 
o 

equal  to  1.0  msec  (t   =  69.3) ,  the  time  constant  used  in 

case  B.   The  nondimensional  damping  coefficient  was  again 

-5 
taken  equal  to  3.416  x  10   ;  and  three  computer  runs  were 

carried  out  at  peak  pressures  of  90  psi,  100  psi,  and  110 
psi.   The  results  are  summarized  in  Table  VII.   Note  that 
for  all  three  runs,  the  integration  was  carried  out  long 
enough  to  include  the  maximum  amplification  of  mode  thir- 
teen. 

Comparison  of  the  amplifications  of  the  first  run  of 
case  B  (Table  V)  with  the  second  run  of  case  D  (Table  VII) 
shows  that  the  amplifications  of  the  Mathieu  modes  were 
again  reduced  approximately  25  to  35  per  cent.  However,  a 
similar  comparison  of  the  third  run  results  of  case  D  with 
the  second  run  results  of  case  B  reveals  that  at  the  higher 
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TABLE  VII 

CASE  D,  AMPLIFICATIONS  AND  STRESSES 

t   =  1.0  msec  (t   =  69.3),  DAMPING  INCLUDED 
o  o 


P   =  90  psi  :  p   =  8.19  x  10~4 

Maximum  circumferential  stress  =  32,100  psi 

Maximum  circumferential  stress  occurred  at  t  =  187 

Integration  terminated  at  T  =  224 

mode      maximum  amp.       t  for  maximum  amp 


221.5 
28.71 
52.59 
58.95 
46.70 
40.03 
26.39 
19.56 
214.7 
158.2 
89.83 
49.03 


0 

1 

26.14 

4 

2.44 

5 

19.15 

6 

61.94 

7 

53.41 

8 

24.94 

9 

11.76 

10 

6.47 

11 

4.66 

13 

199.9 

14 

155.5 

15 

10.75 

P 
o 

=  100  psi 

:  p 

=  9.: 

L0  x 

10" 

■4 

Maximum  circ 

umf erential 

stress 

=  37, 

,900 

psi 

Max 

imum  circ 

:umf  erential 

stress 

occurred 

at  t 

= 

182 

Int 

egration 

terminated  at  t  = 

193 

mode 

maximum 

amp. 

T 

for 

maximum  amp . 

0 

1 

37.87 

4 

3.74 

5 

30.82 

6 

125.6 

7 

113.6 

8 

50.01 

9 

18.98 

10 

9.46 

11 

6.23 

13 

236.5 

14 

208.3 

15 

34.67 

150.2 
162.2 
58.60 
65.44 
52.68 
40.28 
26.93 
19.84 
182.4 
152.5 
83.39 
56.14 
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TABLE  VII   Continued 


P   =  110  psi  :  p   =  lO.Olx  10  4 

o        ^  ^o 

Maximum  circumferential  stress  =  45,400 

Maximum  circumferential  stress  occurred  at  t  =  120 


Integration  terminated  at  t  =  1 


0 

1 

56.98 

4 

7.57 

5 

52.84 

6 

267.9 

7 

245.8 

8 

103.1 

9 

33.80 

10 

13.54 

11 

9.88 

12 

68.58 

13 

257.4 

14 

252.1 

15 

100.5 

mode      maximum  amp.       t  for  maximum  amp 


144.0 
185.9 
63.71 
71.32 
59.19 
46.47 
33.52 
20.08 
176.5 
183.6 
146.5 
118.8 
62.52 


peak  pressure  damping  was  not  nearly  so  effective  in  re- 
ducing the  amplifications  of  the  Mathieu  modes. 

In  Figure  5.13  the  maximum  circumferential  stress  is 
plotted  against  the  peak  pressure.   Assuming  that  the  yield 
point  stress  is  42  ksi,  the  critical  peak  pressure  is 
approximately  106.5  psi.   At  this  pressure  the  amplitudes 
of  both  types  of  modes  are  approximately  the  same.   Thus 
the  Mathieu  modes  play  a  dominant  part  in  the  stress  state 
of  the  cylinder  for  this  loading  condition. 

Critical  Loading  Based  on  Maximum  Stress.   A  critical 
loading  curve  based  on  the  maximum  circumferential  stress 
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at  9  =  0   being  equal  to  the  yield  point  stress  is  plotted 
in  the  pressure-impulse  plane  in  Figure  5.14.   The  curve 
for  an  amplification  of  1000  obtained  by  Anderson  and 
Lindberg  [19]  is  also  shown.   Although  cases  C  and  D 
provide  only  two  points  to  define  the  curve,  clearly, 
critical  stresses  occur  at  pressures  and  impulses  less  than 
those  required  to  produce  amplifications  of  1000. 

Axial  Wave  Form  of  Responding  Modes.   Anderson  and 
Lindberg  [19]  assumed  the  axial  variation  of  the  radial 
displacements  of  the  flexural  modes  and  the  initial  imper- 
fection modes  to  be  a  half  sine  wave.   The  axial  variation 
of  the  displacements  was  completely  unspecified  in  the 
present  analysis.   For  the  solutions  of  this  chapter,  the 
initial  imperfection  modes  were  assumed  to  vary  as  a  half 
sine  wave  in  the  axial  direction,  corresponding  with 
Anderson  and  Lindberg.   In  Figure  5.15,  the  axial  distri- 
butions of  the  computed  radial  displacements  in  modes  6 
and  14  are  plotted  at  arbitrarily  chosen  times.   Note  that 
the  axial  variation  of  mode  6  is  essentially  a  sine  wave 
as  assumed  by  Anderson  and  Lindberg.   On  the  other  hand, 
the  radial  displacement  of  the  Mathieu  mode  is  not  a 
single  half  sine  wave. 

The  axial  variation  of  the  Mathieu  mode  shown  in 
Figure  5.15  appears  to  contain  a  component  of  the  form 
sin(37T^/i)  as  well  as  a  component  of  the  form  sin  (ttC/i  )  • 

The  natural  frequency  Ci        of  a  mode  of  the  form 

^     J      qn 
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wM   sin  (-3— -)  cos  ni 

i 


qn   .  /QTrC 
/^   sin(-a— 

is  given  by 


9        10    9     r,2^2  (1"V2)(^) 

nf-  JL  a^(n^  +  2_JL_)  +  -4__  (5.19) 


qn    12  "  V11       2  '      -   „2  2  2 

(n2+  9LX_) 


A  few  quick  calculations  show  that  when  q  is  equal  to 
three,  the  natural  frequencies  are  equal  to 


n 


ft3n 


10  .3931 

11  .4391 

12  .4964 

13  .5630 

14  .6373 

The  frequency  for  n  =  12  is  very  close  to  one-half,  and 
the  possibility  exists  that  this  mode  might  be  para- 
metrically  excited  by  the  oscillation  of  the  axisymmetric 
mode . 

Random  Imperfections.   In  practice,  the  engineer 
does  not  know  the  magnitude  of  the  initial  imperfection 
until  after  the  shell  has  been  constructed.   This  leads 
one  to  examine  the  sensitivity  of  the  solution  to  the  magni- 
tudes of  the  initial  imperfections.   One  computer  run  was 

n  -5 

made  at  P   =70  psi ,  t   =  3  msec,  and  wi   <  1.0  x  10 
o      ^     o  — 


12 

Mclvor  and  Lovell  [20]  recognized  this  possibility 

and  included  similar  Mathieu  modes  in  their  calculations. 
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The  magnitude  of  the  imperfection  in  each  mode  was  selected 
from  a  table  of  random  numbers.   The  results  obtained  for 
uniform  imperfections  and  random  imperfections  are  shown  in 
Table  VIII  for  purposes  of  comparison. 

Note  that  the  amplifications  of  modes  four  through 
ten  are  nearly  the  same  in  both  cases.   However,  the  dis- 
placements are  significantly  different.   On  the  other  hand, 
the  magnitude  of  the  displacement,  instead  of  the  amplifi- 
cation, remains  approximately  constant  for  modes  13,  14, 
and  15.   Thus,  it  appears  that  the  hyperbolic  modes  are 
very  sensitive  to  the  magnitude  of  the  initial  imperfection, 
whereas  the  response  of  the  Mathieu  modes  is  largely  in- 
dependent of  the  initial  imperfection.   This  is  in  agree- 
ment with  the  conclusions  of  Anderson  and  Lindberg  [19] 
and  Mclvor  and  Lovell  [20].   Anderson  and  Lindberg  con- 
sidered the  hyperbolic  modes  to  be  directly  proportional 
to  the  imperfection;  while  Mclvor  and  Lovell  found  that 
the  maximum  amplitude  of  a  Mathieu  mode  was  rather  in- 
sensitive to  the  initial  imperfection. 

Since  the  stresses  depend  on  the  magnitude  of  the 
displacements,  the  stresses  are  very  sensitive  to  the 
magnitude  of  the  initial  imperfection  when  the  hyperbolic 
modes  play  a  dominant  role. 

Significance  of  Intermodal  Coupling  of  the  Flexural 
Modes .   The  maximum  amplifications  of  the  hyperbolic  modes 
obtained  by  Anderson  and  Lindberg  [19]  agree  well  with  the 
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TABLE  VIII 

EFFECT  OF  MAGNITUDE  OF  INITIAL  IMPERFECTION3 

P   =  70  psi ,  t   =3.0 
o       r      O 

ZERO  DAMPING 


Mode 

.  n 

VI 

Maximum 
Amp . 

i  n.    xlO4 
1 w  1  max 

0 

0.0 

1 

1.0 

x  10"5 

7. 

,93 

.79 

4 

1.0 

x  10~5 

2. 

,22 

.22 

5 

1.0 

x  10~5 

26. 

,70 

2.67 

6 

1.0 

x  10~5 

222. 

,3 

22.23 

7 

1.0 

x  10~5 

101. 

,7 

10.17 

8 

1.0 

x  10~5 

22. 

,02 

2.20 

9 

1.0 

x  10"5 

7. 

,84 

.78 

10 

1.0 

x  10~5 

4. 

,15 

.42 

13 

1.0 

x  10~5 

3. 

,80 

.38 

14 

1.0 

x  10"5 

229. 

,2 

22.92 

15 

1.0 

x  10"5 

4. 

,06 

.41 

Mode 

\: 

.  n 

l 

Maximum 
Amp . 

i  n,    xlO4 
w   max 

0 

0.0 

— 

1 

1.05 

x  10"6 

15. 

,87 

.02 

4 

4.66 

x  10"6 

2. 

,02 

.09 

5 

2.25 

x  10"6 

26. 

,53 

.60 

6 

6.17 

x  10"6 

217. 

,2 

13.40 

7 

6.10 

x  10~7 

102. 

,1 

.623 

8 

5.35 

x  10"6 

22. 

,0 

1.18 

9 

7.12 

x  10"6 

7. 

,83 

.56 

10 

5.69 

x  10"6 

4. 

,15 

.24 

13 

2.53 

x  10"6 

16. 

,70 

.42 

14 

1.80 

x  10"6 

1089. 

,0 

19.6 

15 

1.55 

x  10"6 

3  . 

,57 

.06 

Integration  terminated  at  t  =  197 
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amplifications  computed  by  the  fully  coupled  nonlinear 
theory.   Since  their  analysis  neglects,  among  other  things, 
the  intermodal  coupling  of  the  flexural  modes,  the  possi- 
bility arises  that  the  coupling  between  the  flexural 
modes  is  of  negligible  importance.   Neglecting  intermodal 
coupling  of  the  flexural  modes  is  equivalent  to  retaining 
only  those  nonlinear  products  which  contain  the  Fourier 
coefficient  of  the  axisymmetric  mode.   Motion  in  the 
axisymmetric  mode  affects  the  motion  of  each  of  the  other 
modes,  and  conversely  each  mode  couples  back  to  the  axi- 
symmetric mode.   However,  modes  i  and  j  do  not  affect  each 
other  if  neither  i  nor  j  is  equal  to  zero. 

The  computer  program  can  be  altered  very  readily  to 
effect  this  limited  coupling.   A  run  was  made  for  the 
parameters  of  case  A,  P   =70  psi,  using  the  altered 
program.   The  maximum  amplifications  for  both  complete 
and  limited  nonlinear  coupling  are  presented  in  Table  IX. 
The  integration  was  terminated  at  r  =  200  so  that  mode 
13  did  not  have  time  to  grow.   The  magnitudes  and  times 
of  the  maximum  amplification  of  all  of  the  modes  are  nearly 
the  same  for  both  cases.   On  the  basis  of  this  one  solu- 
tion, it  appears  that  the  intermodal  coupling  of  the 
flexural  modes  for  this  problem  is  quite  weak  and  can  be 
neglected . 

Additional  computer  runs  should  be  made  in  the  future 
to  cover  a  larger  spectrum  of  loading  conditons.  If  it  can 
be  shown  that  the  intermodal  coupling  of  the  flexural  modes 
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p  = 


TABLE  IX 

AMPLIFICATIONS  WITH  TOTAL 
AND  LIMITED  INTERMODAL  COUPLING 
70  psi,  t 


o 


n  -5 

3.0  msec,  wi   =  1.0  x  10 


ZERO  DAMPING 


Total  Coupl ing 

Limited  Coupling 

Mode 

i 
Maximum 

Amplification 

t  For 
Max.  Amp. 

Maximum 
Amplification 

t  For 
Max.  Amp. 

0 

1 

4 

5 

6 

7 

8 

9 

10 

13 

14 

15 

7.933 
2.220 
26.70 

222.3 

101.7 
22.02 
7.  838 
4.147 
3.796 

229.2 
4.064 

190.6 
88.35 
71.71 
97.35 
72.23 
40.37 
26.41 
19.43 
193.8 
164.6 
162.8 

.039 
1.974 
26.20 
219.6 
101.3 
21.98 
7.827 
4.143 
8.854 
229.  3 
3.795 

198.1 
33.53 
71.78 
97.48 
72.15 
40.37 
26.56 
19.45 
199.9 
164.6 
191.4 

Max. 
Stress 

30,000  psi 

@  t  =  19  7 

32,900  psi 

@  x  =  165 

165 


is  indeed  negligible  for  all  loadings  of  interest,  the 
analysis  and  the  computer  program  can  be  simplified  very 
substantially . 

Significance  of  Coupling  Terms  in  Mode  Zero  Equations 
The  coupling  terms  in  the  equations  of  mode  0,  the  axi- 
symmetric  mode,  provide  the  mechanism  for  the  extraction 
of  energy  from  that  mode.   A  natural  question  to  ask  is 
whether  the  response  of  the  flexural  modes  would  be 
affected  if  this  mechanism  were  removed.   If  this  coupling, 
as  well  as  the  intermodal  coupling  of  the  flexural  modes, 
is  unimportant,  then  the  response  of  all  of  the  flexural 
modes  could  be  computed  one  mode  at  a  time,  using  equations 
(5.5)  and  (5.6),  or  similar  equations.   Three  computer 
runs  were  made  to  study  this  question.   In  all  three  runs, 
the   coupling  terms  in  the  mode  0  equations  were  equated 
to  zero.   Run  one  included  the  flexural  modes  5  through  8, 
run  two  included  the  flexural  modes  12  through  15,  and  run 
three  included  the  flexural  modes  16  through  18.   The 
choice  of  these  flexural  modes  rules  out  any  intermodal 
coupling  of  the  flexural  modes.   The  loading  had  a  peak 
pressure  of  70  psi  and  a  time  constant  of  3.0  msec,  the 
same  as  the  loading  of  the  second  run  of  case  A.   The 
maximum  amplification  of  each  mode  from  both  solutions  is 
contained  in  Table  X.   Comparison  of  the  appropriate 
amplifications  from  Table  IV  with  those  contained  in  Table 
X  reveals  that  the  amplifications  of  the  hyperbolic  modes, 
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TABLE  X 


RESPONSE  OF  INDIVIDUAL  MODES  WITHOUT 

ENERGY  EXTRACTION  FROM  THE  AXISYMMETRIC  MODE 

P   =  70  psi  (p   =  6.37  x  10~4),  t   =  3.0  msec  (x   =  208) 
o       c  ^o  o  o 


.„  ,  Maximum  Amplification 

Run         Mode  _.       ,  c .     „         -.  .  ~ 

Attained  at  T  =  140 


5  26.8 

6  228.0 

7  104.0 

8  21.4 

12  1.8 

13  1.8 

14  1990.0 

15  3.4 

16  1.7 

17  1.7 

18  2.7 
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modes  5  through  8,  were  not  changed  appreciably  when  the 
coupling  terms  in  the  mode  0  equation  were  neglected.   On 
the  other  hand,  the  amplification  of  the  Mathieu  mode, 
mode  14,  was  drastically  larger.   This  is  reasonable  since 
the  initial  conditons  rather  than  the  loading  cause  the 
oscillations  of  the  axisymmetric  mode.   Thus,  the  oscilla- 
tory energy  of  the  axisymmetric  mode  is  limited.   If  this 
limited  energy  can  not  be  removed  by  coupling,  then  the 
Mathieu  mode  will  grow  very  large  as  if  an  oscillatory 
load  were  being  applied. 

Although  the  results  are  somewhat  limited  in  number, 
they  indicate  that  the  response  of  a  particular  flexural 
mode  has  little  direct  effect  on  the  response  of  another 
flexural  mode;  that  is  the  intermodal  coupling  of  the  flex- 
ural modes  is  weak.   However,  all  of  the  flexural  modes 
have  a  very  significant  indirect  influence  on  the  response 
of  the  Mathieu  modes,  since  they  extract  energy  from  the 
axisymmetric  mode.   Consequently  the  response  of  all  the 
modes  must  be  computed  simultaneously  in  order  to  obtain 
the  total  response  of  the  cylindrical  shell. 
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CHAPTER  VI 
SUMMARY,  FINDINGS,  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

I .   SUMMARY 

A  numerical  algorithm  was  developed  for  computing 
the  dynamic  response  of  a  ring-stiffened,  nearly  circular 
cylindrical  shell  under  transient  axisymmetric  loads. 
The  governing  equations,  which  are  based  on  the  nonlinear 
Donnell  theory,  account  for  the  class  of  radial  imperfec- 
tions whose  circumferential  variation  can  be  expressed 
in  a  Fourier  cosine  series.   The  equations  for  the  rings 
and  shell  structure  were  derived  by  using  Hamilton's 
principle  and  include  the  radial  inertia  and  a  viscous 
damping  term.   The  analysis  accounts  for  the  discrete 
location  and  eccentricity  of  a  set  of  identical  rings. 

Elimination  of  the  circumferential  variable  was 
achieved  by  assuming  a  Fourier  sine  or  cosine  series  for 
the  displacements.   The  nonlinear  terms  generated  products 
of  Fourier  series  which  were  expanded  to  single  series, 
yielding  modal  coupling  terms.   The  series  were  truncated 
for  the  actual  computations,  with  a  maximum  of  fifteen 
participating  modes.   The  retained  modes  were  not  limited 
to  the  first  fifteen,  but  were  selected  to  include  those 
modes  most  responsive  to  the  particular  load.   The  partial 
differential  equations  were  replaced  with  difference- 
differential  equations  by  approximating  the  axial 
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derivatives  with  finite  differences.   A  nonconventional 
differencing  procedure,  commonly  called  "modified"  or 
"half  station"  differencing,  was  used.   Provisions  were 
made  for  both  the  clamped  and  simply  supported  end 
conditions.   The  iterative  Newmark  beta-method  was  employed 
for  the  time  integration  of  the  radial  difference- 
differential  equations.   The  axial  and  circumferential 
equilibrium  equations  were  solved  by  Gauss  elimination 
during  each  iteration. 

The  algorithm  was  used  to  compute  the  response  of 
a  simply  supported  shell  under  a  uniform,  exponentially 
decaying,  radial  pressure.   This  particular  shell  and  simi- 
lar loading  conditions  have  been  considered  in  the  earlier 
investigations  of  Anderson  and  Lindberg  [19]  and  Mclvor 
and  Lovell  [20].   The  load  parameters  were  selected  so  as 
to  fall  in  the  regime  where  dynamic  buckling  is  governed 
by  an  elastic  model,  and  the  peak  pressures  and  total 
impulses  considered  were  between  the  static  buckling  limit 
and  the  dynamic  buckling  limit.   Several  response  histories 
were  computed,  both  with  and  without  viscous  damping;  and 
comparisons  were  made  with  the  results  of  References 
[19,  20].   The  significance  of  the  intermodal  coupling  of 
the  flexural  modes  and  the  feedback  coupling  to  the  axi- 
symmetric  mode  was  also  investigated.   Finally,  the  role 
played  by  the  imperfection  magnitudes  was  briefly  considered 
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A  potential  user  of  the  computer  program,  or  one 
interested  in  pursuing  this  work  further ,  should  be  aware 
of  the  following  limitations  of  the  program: 

1.  The  analysis  employs  the  Donnell  shell  theory 

which  is  not  valid  for  many  classes  of  shells 
[48]. 

2.  The  circumferential  variation  of  the  initial 

imperfections  is  expressed  in  a  Fourier  cosine 
series,  instead  of  a  complete  series  containing 
both  sines  and  cosines. 

3.  The  axial,  circumferential,  and  rotatory  inertia 

of  the  rings  and  the  shell  are  neglected. 

4.  The  portions  of  the  computer  program  that  cal- 

culate ring  effects  have  been  verified  in  only 
a  cursory  manner. 

5 .  A  maximum  of  three  rings  can  be  accomodated  by 

the  computer  program.  The  rings  must  be  iden- 
tical in  geometry  and  material  properties  and 
spaced  so  as  to  divide  the  shell  into  segments 
of  equal  length.  A  principal  axis  of  the  ring 
cross  sections  must  be  parallel  to  the  axis  of 
the  cylinder. 

6.  The  imperfection  of  the  shell  must  not  distort 

the  rings. 

7.  The  rotations  of  the  rings  must  be  small,  and 

restrictive  warping  is  not  accounted  for  in  the 
expression  for  the  strain  energy  of  the  rings. 
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8.  The  distance  between  the  midsurface  of  the  shell 

and  the  centroids  of  the  ring  cross  sections 
must  be  small  compared  to  the  radius  of  the 
cylinder. 

9.  The  circumferential  stress  is  computed  only  at 

6=0.   The  meridional  curvature  change  is 

not  included  in  the  calculation  of  the  cir- 

.  .  2 

cumferential  stress;  in  addition,  1  -v    is 

taken  equal  to  unity. 

10.   In  the  axial  equilibrium  equation,  the  effect 

of  a  ring  is  felt  at  a  station  which  is 

located  one-half  of  the  mesh  spacing  from  the 

ring. 

II.   FINDINGS 

The  exponentially  decaying  pressure  pulses  were  found 
to  excite  flexural   modes  belonging  to  two  families,  herein 
referred  to  as  hyperbolic  and  Mathieu  modes.   The  hyperbolic 
modes  are  the  buckling  modes,  and  they  grow  to  their  maximum 
amplitude  in  an  exponential  manner  with  little  accompanying 
oscillation.   On  the  other  hand,  the  Mathieu  modes  grow  in 
an  oscillatory  manner  requiring  several  cycles  to  attain 
their  maximum  amplitude. 

The  hyperbolic  modes  grow  primarily  as  a  result  of 
the  quasi-static  displacement  of  the  axisymmetric  mode. 
The  loads  considered  in  the  present  study  were  on  the  shell 
for  a  period  of  time  much  longer  than  the  natural  period 
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of  the  axisymmetric  mode.   Hence,  the  quasi-static  response 
of  the  axisymmetric  mode  was  nearly  proportional  to  the 
instantaneous  magnitude  of  the  pressure.   A  given  hyper- 
bolic mode  was  found  to  respond  as  long  as  the  instan- 
taneous magnitude  of  the  pressure  exceeded  the  static 
buckling  pressure  of  the  mode. 

The  Mathieu  modes  were  found  to  be  those  modes  whose 
natural  frequencies  were  approximately  one-half  of  the 
axisymmetric  natural  frequency  and  usually  much  higher 
than  the  frequencies  of  the  hyperbolic  modes.   The  Mathieu 
modes  were  identified  through  the  use  of  a  modified 
Mathieu  stability  chart.   The  modification  replaced  the 
natural  frequency  of  the  flexural  mode  by  a  pseudo  fre- 
quency which  accounted  for  the  instantaneous  value  of 
the  pressure.   Hence  the  susceptibility  of  a  given  flex- 
ural mode  to  parametric  excitation  depended  not  only  on 
the  peak  pressure  and  frequencies  of  the  flexural  and 
axisymmetric  modes,  but  also  on  the  instantaneous  value 
of  the  pressure.   Thus,  as  the  pressure  pulse  decayed,  the 
points  corresponding  to  the  flexural  modes  moved  across 
the  Mathieu  diagram.   Some  modes  moved  from  an  unstable 
region  into  a  stable  region,  while  other  modes  moved  from 
a  stable  region  into  an  unstable  region.   The  response  of 
the  Mathieu  modes  could  be  characterized  as  beat-type 
oscillations  during  which  energy  was  cyclically  exchanged 
between  modes. 
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Anderson  and  Lindberg  [19]  investigated  the  dynamic 
pulse  buckling  of  cylindrical  shells  subjected  to  axi- 
symmetric  loadings.   They  showed  that  dynamic  buckling 
occurs  as  a  result  of  the  rapid  growth  of  the  hyperbolic 
modes,  and  that  the  Mathieu  modes  are  of  little  importance. 

Mclvor  and  Lovell  [20]  studied  the  nonlinear  dynamic 
response  of  shells  under  purely  impulsive  loads.   Their 
results  show  that  for  such  loads  only  the  Mathieu  modes 
respond  since  the  quasi-static  response  of  the  axisymmetric 
mode  is  zero.   However,  the  stresses  developed  during  the 
response  were  found  to  be  several  times  larger  than 
stresses  of  an  axisymmetric  response  to  the  same  impulsive 
load . 

The  results  of  the  present  work  consolidate  the  find- 
ings of  these  two  previous  investigations.   At  the  lower 
peak  pressures  and  total  impulses,  for  which  the  quasi- 
static  response  of  the  axisymmetric  mode  was  small,  the 
Mathieu  modes  dominated;  and  the  computed  motion  was  quali- 
tatively similar  to  the  motion  computed  by  Mclvor  and 
Lovell.   As  the  peak  pressure  and  total  impulse  were 
increased,  a  gradual  transition  occurred;  and  the  hyper- 
bolic modes  played  an  increasingly  significant  role.   Near 
the  dynamic  buckling  threshold,  the  hyperbolic  modes  com- 
pletely dominated  the  Mathieu  modes.   The  dynamic  buckling 
loads  predicted  by  the  fully  coupled  nonlinear  analysis 
agreed  very  closely  with  Anderson  and  Lindberg 's  results. 
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However,  at  peak  pressures  and  total  impulses  well  below 
the  dynamic  buckling  threshold,  stresses  in  excess  of  the 
yield  point  stress  occurred  due  to  the  combined  response 
of  the  Mathieu  and  hyperbolic  modes. 

For  those  solutions  in  which  a  small  amount  of  damping 
was  assumed,  the  maximum  amplitudes  of  both  the  dominant 
hyperbolic  and  Mathieu  modes  were  reduced  by  10  to  30  per 
cent  in  comparison  to  the  zero  damping  case.   Generally, 
the  Mathieu  modes  were  attenuated  more  than  the  hyperbolic 
modes.   A  critical  pressure-impulse  curve  based  on  the 
yield  point  stress  was  determined  for  the  damped  response. 
This  curve  shows  that  stresses  in  excess  of  the  yield  point 
stress  can  occur  at  total  impulses  and  peak  pressures  less 
than  those  required  for  dynamic  buckling  even  when  damping 
is  accounted  for. 

The  intermodal  coupling  of  the  flexural  modes  and  the 
feedback  coupling  to  the  axisymmetric  mode  did  not  signif- 
icantly affect  the  response  of  the  hyperbolic  modes. 
Furthermore,  the  maximum  amplitude  attained  by  a  given 
hyperbolic  mode  was  found  to  be  linearly  proportional  to 
shape  and  magnitude  of  the  imperfection  in  that  mode.   These 
findings  are  in  complete  agreement  with  the  assumptions  made 
by  Anderson  and  Lindberg  [19]. 

However,  the  feedback  coupling  to  the  axisymmetric 
mode  very  significantly  affected  the  response  of  the  Mathieu 
modes.   This  coupling  is  the  mechanism  by  which  energy  is 
exchanged  between  the  axisymmetric  mode  and  the  flexural 
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modes.   On  the  other  hand,  the  maximum  amplitudes  attained 
by  the  Mathieu  modes  were  unaffected  by  both  the  inter- 
modal  coupling  of  the  flexural  modes  and  the  magnitude  of 
the  imperfections . 

The  findings  presented  here  apply  to  one  specific 
simply  supported  shell  and  loading  condition.   The  buckling 

TT  X 

mode  shape  of  this  shell  is  sin^-cos66,  and  the  dominant 

J_i 

Mathieu  mode  is  sin^p-cosl40 .  For  other  shells  with  differ- 
ent  loading  conditons,  hyperbolic  modes,  and  Mathieu  modes, 
the  coupling  may  be  more  significant;  and  these  conclusions 
would  not  apply.   This  remains  to  be  seen. 

The  algorithm  developed  in  this  dissertation  is 
capable  of  computing  the  fully  coupled,  nonlinear  dynamic 
response  of  cylindrical  shells  under  transient  pressures. 
However,  the  practicality  of  the  method  is  somewhat  limited 
by  the  execution  time  of  the  computer  program.   For  example, 
approximately  forty-five  minutes  are  required  on  the 
IBM  360  model  67  to  compute  a  solution  through  thirty 
cycles  of  the  axisymmetric  mode  when  fifteen  Fourier  modes 
are  retained.   When  the  pressure  pulse  decays  rather  slowly, 
as  it  did  in  cases  B  and  D  of  Chapter  V,  the  solution  must 
be  computed  through  as  many  as  fifty  cycles  of  the  axi- 
symmetric mode  in  order  to  encompass  the  response  of  all 
of  the  Mathieu  modes. 
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III.   RECOMMENDATIONS  FOR  FUTURE  WORK 

The  computer  program  has  yet  to  be  fully  exploited; 
and,  as  previously  pointed  out,  the  correctness  of  those 
portions  of  the  program  which  compute  the  ring  contributions 
have  not  been  fully  confirmed.   Consequently,  the  following 
studies  should  be  conducted  using  the  program: 

1.  Additional  solutions  carried  out  to  longer  times 

are  needed  to  establish  the  significance  of 
the  intermodal  coupling  to  the  response  of  the 
Mathieu  modes.   If  part  of  the  coupling  terms 
could  be  eliminated,  the  computing  time  could 
be  reduced  appreciably,  thus,  extending  the 
usefulness  of  the  program. 

2.  Solutions  should  be  generated  for  shorter  shells 

with  L/a  =  1.   For  shells  of  this  geometry, 
the  static  buckling  mode  and  the  Mathieu  mode 
coincide.   The  dynamic  buckling  loads  might 
turn  out  to  be  significantly  different  from 
those  predicted  by  Reference  [19] . 

3.  Solutions  should  be  computed  for  shells  with 

clamped  ends  and  compared  with  the  results 
for  the  simply  supported  shell. 

4.  The  stabilizing  or  destabilizing  effects  of  ring 

stiffeners  should  be  studied.   For  example,  if 
the  stiffeners  are  too  close  together,  it  is 
possible  that  coincidence  of  the  Mathieu  modes 
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and  the  hyperbolic  modes  may  render  the  struc- 
ture unstable. 

The  program  should  be  used  to  study  the  dynamic 
stability  of  shells  under  moving  loads. 

An  effort  should  be  made  to  improve  the  Fortran 
coding  of  the  program.   It  is  believed  that  a 
professional  programmer  could  improve  the 
coding  and,  thus,  reduce  the  running  time. 


178 


LIST  OF  REFERENCES 


1.  Culbertson,  D.  W. ,  "Conical  Shock  Tube  Blast  Simulator 
Facility  Description  and  Capabilities,"  NWL-TM  No. 
T-15/67,  U.S.  Naval  Weapons  Laboratory,  Dahlgren, 
Virginia. 

2.  Bolotin,  V.  V.,  The  Dynamic  Stability  of  Elastic 
Systems ,  translation  from  the  Russian  (1956)  ,  Holden- 
Day,  Inc.,  San  Francisco,  London,  Amsterdam.  1964. 

3.  Agamirov,  V.  L.  and  Vol'mir,  A.  S.,  "Behavior  of 
Cylindrical  Shells  Under  Dynamic  Loading  by  Hydrostatic 
Pressure  or  by  Axial  Compression,"  translation  from 
the  Russian  (1959),  ARS  Journal,  Vol.  31,  No.  1,  Jan. 
1961,  pp.  98-102. 

4.  Yao ,  J.  C. ,  "Stability  of  a  Cylinder  Under  Dynamic 
Radial  Pressure,"  ARS  Journal,  Vol.  31,  No.  12,  Dec. 
1961,  pp.  1705-1708. 

5.  Kadashevich,  Y.  I.  and  Pertsev,  A.  K. ,  "Loss  of  Stabil- 
ity of  a  Cylindrical  Shell  Under  Dynamic  Loads,"  trans- 
lation from  the  Russian  (1960),  ARS  Journal ,  Vol.  32, 
No.  1,  Jan.  1962,  pp.  140-143. 

6.  Yao,  J.  C. ,  "Dynamic  Stability  of  Cylindrical  Shells 
Under  Static  and  Periodic  Axial  and  Radial  Loads," 
AIAA  Journal,  Vol.  1,  No.  6,  Jun.  1963,  pp.  1391-1396. 

7.  McLachlan,  N.  W. ,  Theory  and  Application  of  Mathieu 
Functions ,  Oxford  University  Press,  New  York,  London, 
1947. 

8.  Wood,  J.  D.  and  Koval ,  L.  R. ,  "Buckling  of  Cylindrical 
Shells  under  Dynamic  Loads,"  AIAA  Journal ,  Vol.  1, 

No.  11,  Nov.  1963,  pp.  2576-2582. 

9.  Goodier,  J.  N.  and  Mclvor,  I.  K. ,  "The  Elastic  Cylin- 
drical Shell  Under  Nearly  Uniform  Radial  Impulse," 
Journal  of  Applied  Mechanics,  Transactions  of  the 
ASME,  Vol.  31,  No.  2,  Jun.  1964,  pp.  259-266. 

10.   Lindberg,  H.  E. ,  "Buckling  of  a  Very  Thin  Cylindrical 

Shell  Due  to  an  Impulsive  Pressure,"  Journal  of  Applied 
Mechanics ,  Transactions  of  the  ASME,  Vol.  31,  No.  2, 
Jun.  1964,  pp.  267-272. 


179 


11.  Roth,  R.  S.  and  Klosner,  J.  M. ,  "Nonlinear  Response 
of  Cylindrical  Shells  Subjected  to  Dynamic  Axial 
Loads/'  AIAA  Journal,  Vol.  2,  No.  10,  Oct.  1964, 
pp.  1788-1794. 

12.  Evan-Iwanowski ,  R.  M. ,  "On  the  Parametric  Response 
of  Structures,"  Applied  Mechanics  Reviews,  Vol.  18, 
No.  9,  Sept.  1965,  pp.  699-702. 

13.  Bienick,  M.  P.,  Fan,  T.  C. ,  and  Lackman,  L.  M. , 
"Dynamic  Stability  of  Cylindrical  Shells,"  AIAA 
Journal ,  Vol.  4,  No.  3,  Mar.  1966,  pp.  495-500. 

14.  Hegemeir,  G.  A.,  "Instability  of  Cylindrical  Shells 
Subjected  to  Axisymmetric  Moving  Loads,"  Journal  of 
Applied  Mechanics,  Transactions  of  the  ASME,  Vol.  33, 
No.  2,  Jun.  1966,  pp.  289-296. 

15.  Hegemeir,  G.  A.,  "Stability  of  Cylindrical  Shells 
Under  Moving  Loads  by  the  Direct  Method  of  Liapunov," 
Journal  of  Applied  Mechanics,  Transactions  of  the 
ASME,  Vol.  34,  No.  4,  Dec.  1967,  pp.  991-998. 

16.  Fan,  T.  C,  "Dynamic  Stability  of  a  Thin  Cylinder 
Under  Radial  Pressure,"  Memorandum  RM-4583-PR,  Jun. 
1966,  The  Rand  Corp. 

17.  Lock,  M.  H. ,  "The  Nonlinear  Response  of  a  Thin 
Circular  Ring  Under  a  Uniform  Step  Pressure  Load," 
TR-1001 (2240-30)-l,  Nov.  1966,  Aerospace  Corp. 

18.  Lindberg,  H.  E.,  Anderson,  D.  L.,  Firth,  R.  D. ,  and 
Parker,  L.  V.,  "Response  of  Reentry  Vehicle  -  Type 
Shells  to  Blast  Loads,"  Final  Report  to  Lockheed 
Missiles  and  Space  Co.,  Contract  P.O.  24-14517  under 
AF  04(694)  -  655,  1965,  Stanford  Research  Institute. 

19.  Anderson,  D.  L.  and  Lindberg,  H.  E.  "Dynamic  Pulse 
Buckling  of  Cylindrical  Shells  Under  Transient  Lateral 
Pressures,"  AIAA  Journal ,  Vol.  6,  No.  4,  Apr.  1968, 
pp.  589-598. 

20.  Mclvor,  I.  K.  and  Lovell,  E.  G.,  "Dynamic  Response  of 
Finite-Length  Cylindrical  Shells  to  Nearly  Uniform 
Radical  Impulse,"  AIAA  Journal,  Vol.  6,  No.  12,  Dec. 
1968,  pp.  2346-2351. 

21.  Dietz,  W.  K. ,  "On  the  Dynamic  Stability  of  Eccen- 
trically Reinforced  Circular  Cylindrical  Shells," 
Ph.D.  dissertation,  Jan.  1967,  Syracuse  University. 


180 


22.  Klosner,  J.  M.  and  Franklin,  H.  N.,  "A  Note  on  the 
Dynamic  Instability  of  Stiffened  Cylindrical  Shells," 
PIBAL  Report  No.  68-1,  Jan.  1968,  Polytechnic  Insti- 
tute of  Brooklyn,  Department  of  Aerospace  Engineering 
and  Applied  Mechanics. 

23.  Block,  D.  L. ,  "Influence  of  Discrete  Ring  Stiffeners 
and  Prebuckling  Deformations  on  the  Buckling  of 
Eccentrically  Stiffened  Orthotropic  Cylinders," 
NASA  TN  D-4283,  Jan.  1968. 

24.  Donnell,  L.  H.,  "Stability  of  Thin-walled  Tubes  Under 
Torsion,"  NACA  Report  No.  4  79,  193  3. 

25.  Donnell,  L.  H.,  "A  New  Theory  for  the  Buckling  of  Thin 
Cylinders  Under  Axial  Compression  and  Bending," 
Transactions  of  the  American  Society  of  Mechanical 
Engineers,  Vol.  56,  1934,  pp.  795-806. 

26.  Galletly,  G.  D. ,  "On  the  In-vacuo  Vibrations  of  Simply 
Supported,  Ring-Stiffened  Cylindrical  Shells,"  Pro- 
ceedings of  the  U.S.  National  Congress  of  Applied 
Mechanics,  2nd. ,  1954,  pp.  225-231. 

27.  Goldstein,  H. ,  Classical  Mechanics,  Addison-Wesley , 
Reading,  Massachusetts,  1950,  pp.  38-44. 

28.  Timoshenko,  S.  and  Goodier,  J.  N.,  Theory  of  Elas- 
ticity ,  2nd  ed. ,  McGraw-Hill,  New  York,  Toronto, 
London,  1951,  pp.  275-280. 

29.  Hurty,  W.  C.  and  Rubinstein,  M.  F.,  Dynamics  of  Struc- 
tures ,  Prentice  Hall,  Englewood  Cliffs,  New  Jersey, 
1964,  pp.  80-90. 

30.  Ball,  R.  E.,  "A  Geometrically  Nonlinear  Analysis  of 
Arbitrarily  Loaded  Shells  of  Revolution,"  NASA  CR-909, 
Jan.  1968. 

31.  Chuang,  K.  P.  and  Veletsos ,  A.  S.,  "A  Study  of  Two 
Approximate  Methods  of  Analyzing  Cylindrical  Shell 
Roofs,"  Civil  Engineering  Studies,  Structural  Research 
Series,  No.  258,  University  of  Illinois,  Oct.  1962. 

32.  Ball,  R.  E.,  "Dynamic  Analysis  of  Rings  by  Finite 
Differences,"  Journal  of  the  Engineering  Mechanics 
Division,  Proceedings  of  the  ASCE ,  Vol.  93,  No.  EMI , 
Mar.  22,  1966,  pp.  1-10. 

33.  Cyrus,  N.  J.  and  Fulton,  R.  E.,  "Finite  Difference 
Accuracy  in  Structural  Analysis,"  Journal  of  the 
Structural  Division,  Proceedings  of  the  ASCE,  Vol. 
92,  No.  ST6,  Dec.  1966,  pp.  459-471. 


181 


34.  Wah,  T.  and  Lu,  W.  C.  L.  ,  "Vibration  Analysis  of 
Stiffened  Cylinders  including  Inter-Ring  Motion," 
The  Journal  of  the  Acoustical  Society  of  America, 
Vol.  43,  No.  5,  1968,  pp.  1005-1016. 

35.  Smith,  G.  D. ,  Numerical  Solutions  of  Partial  Differ- 
ential Equations,  Oxford  University  Press,  New  York, 
London,  1965,  pp.  143-161. 

36.  Forsy the ,  G.  E.  and  Wasow,  W.  R. ,  Finite-Difference 
Methods  for  Partial  Differential  Equations,  John 
Wiley  &  Sons,  New  York,  London,  Sydney,  1960, 

pp.  208-220. 

37.  Budiansky,  B.  and  Radkowski,  P.  P.,  "Numerical  Analysis 
of  Unsymmetrical  Bending  of  Shells  of  Revolution," 
AIAA  Journal,  Vol.  1,  No.  8,  Aug.  1963,  pp.  1833-1842. 

38.  Jennings,  W. ,  First  Course  in  Numerical  Methods, 
Macmillan,  New  York,  1964,  p.  149. 

39.  Crandall,  S.  H. ,  Engineering  Analysis,  McGraw-Hill, 
New  York,  Toronto,  London,  1956,  pp.  376-396. 

40.  Johnson,  D.  E.  and  Grief,  R. ,  "Dynamic  Response  of 
a  Cylindrical  Shell:  Two  Numerical  Methods,"  AIAA 
Journal,  Vol.  4,  No.  3,  Mar.  1966,  pp.  486-494. 

41.  Forsythe  and  Wasow,  op.  cit. ,  pp.  101-107. 

42.  Newmark,  N.  M. ,  "Computation  of  Dynamic  Structural 
Response  in  the  Range  Approaching  Failure,"  Proceedings 
of  the  Symposium  on  Earthquake  and  Blast  Effects  on 
Structures,  Earthquake  Engineering  Institute,  Los 
Angeles,  1952. 

43.  Newmark,  N.  M. ,  "A  Method  of  Computation  for  Struc- 
tural Dynamics,"  Journal  of  the  Engineering  Mechanics 
Division,  Proceedings  of  the  ASCE,  Vol"!  85,  No.  EM  3 , 
Jul.  1959,  pp.  67-94. 

44.  Chang,  G.  C.  and  Ang ,  A.  H.  S.,  "Numerical  Analysis  of 
Plane  Structure-Medium  Interaction  in  Elastic  Per- 
fectly Plastic  Media,"  Civil  Engineering  Studies, 
Structural  Research  Series  No.  307,  University  of 
Illinois,  Jun .  1966. 

45.  Bhuta,  P.  G. ,  "Transient  Response  of  a  Thin  Elastic 
Cylindrical  Shell  to  a  Moving  Shock  Wave,"  The  Journal 
of  the  Acoustical  Society  of  America,  Vol.  35,  No.  1, 
Jan.  1963,  pp.  25-30. 


182 


46.  Fung,  Y.  C.  ,  Sechler,  E.  E. ,  and  Kaplan,  A.,  "On  the 
Vibration  of  Thin  Cylindrical  Shells  under  Internal 
Pressure,"  Journal  of  the  Aeronautical  Sciences , 
Vol.  24,  No.  9,  Sep.  1957,  pp.  650-660. 

47.  Stoker,  J.  J.,  Nonlinear  Vibratons ,  Interscience 
Publishers,  Inc.,  New  York,  1950,  pp.  202-213. 

48.  Krauss,  H. ,  Thin  Elastic  Shells,  John  Wiley  &  Sons, 
Inc.,  New  York,  London,  Sydney,  1967,  pp.  221-229. 


183 


APPENDIX  A 
INCONSISTENCY  OF  THE  CONVENTIONAL  DIFFERENCES 

We  will  first  show  that  the  difference  equations 
obtained  by  discretization  of  the  governing  differential 
equations  do  not  agree  with  the  difference  equations  ob- 
tained by  minimizing  the  finite  difference  form  of  the 
potential  energy  when  conventional  central  differences  are 
used  for  approximating  the  derivatives.   We  will  then  show 
that  the  two  sets  of  difference  equations  are  identical 
when  modified  differences  are  used.   Chuang  and  Veletsos 
[31]  presented  a  very  similar  proof  of  this  inconsistency 
of  the  conventional  differences,  but  they  did  not  expand 
the  displacements  in  Fourier  series  as  has  been  done  here. 

Consider  the  static  equilibrium  of  a  cylindrical 

shell  under  a  uniform  radial  pressure.   For  the  purpose 

of  the  present  discussion,  the  linear  shell  equations  will 

be  used.   After  deleting  the  acceleration  and  velocity 

terms,  equations  (3.15)  become 

n      1  +  v     n    1-v   2  n      n     A 
u,c-c-  +  — 2~"  nV'r  ~  ~2~   n  u   -  vw,^  =  0 

2  n    1+v    n     1-v  n        n    n 
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(A.l) 


~         n     n      n 
o  n   p  +  nv   -  w   +  vu ,  r 
0n^  '  t, 
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The  conventional  central  differences  for  the  first,  second, 
and  fourth  derivatives  of  a  function  f (?)  are 

{f'^K   =  2d  (fK+l  "  fK-l) 

(f'^K  =  72  (fK+l  "  2fK  +  fK-l»  (A'2> 

d 

(f'e««'K  =  jt    <fK+2  "  4fK+l  +  6fK  -  4fK-l  +  fK-2» 

where  d  is  spacing  of  the  node  points. 

Using  the  conventional  central  differences  given  by 
(A. 2)  to  approximate  the  derivatives  in  equations  (A.l) 
leads  to 

1  ,  n      „  n     n   >    1+v    1  ,  n     n   .    1-v   2  n 
^2  (UK+1  "  2UK  +  UK-1>  +  —  n  2d(VK+l-VK-l}  "  ~   n  UK 

1   ,  n       n   >    n 
"VH  (WK+1  "  WK-1}  =  ° 

2  n    1+v    1   .  n     n   , ,  1-v  1  .  n    _  n,  n. ,    n   _ 
"n  VK  "  I"  n  ii  (UK+1-UK-1)+  —  3-(VK+l-2VK+VK)+  nWK=  ° 


1   2  r  1   /  n    >•  n  ^n    .  n      n.               ,,  ,, 

12  a  [7  (WK+2"4WK+1+  6WK  "  4WK-1  +  WK-2}  (A'3) 

~  2  1   ,  n      „  n  n   ,  ,   4  n, 

_2n   -2  (wR+   -  2wK  +  v        )    +  n  w  ]  = 

d 

,         _       n     n  ,  1   .  n       n   x 

60n  PK  +  nVK  "  WK  +  V2d  (UK+1  "  UK-1} 


A  set  of  difference  equations  can  also  be  derived  by 
minimizing  the  potential  energy.   The  potential  energy  II, 
for  the  linear  shell  under  consideration  here,  is  given  by 
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1       L      27ia 

200  xx  xy'xy  yy 


-My   +2M   x    "MY   -  2PW)  dydx 
xAx      xyAxy     yAy  ■* 


(A. 4) 


The  strain-displacement  relations  of  the  linear  theory  are 


£       =    U,  X  FT 

xx  A   =  W, 

x      xx 

e   =  V,   -  -  x   =  W,  (A. 5) 

y     y   a         Ay     yy 


Y    =  U,   +  V,         x     =  w' 
1  xy     y     x       Axy      xy 


The  in-plane  forces  and  moments  are  related  to  the 
displacements  by 

W 
N    =  B[U,   +  v(V,   -  -)  ] 
x        'x       y    a 

W 
N    =  B  V,   -  -  +  vU,  ) 
Y         y    a       x 

N    =  (^)B(U,   +  V,  ) 
xy      2       y      x 


(A. 6 


M    =  -D(W,    +  vW,   ) 
x  xx      'yy 


M    =  -D(W,    +  vW,   ) 

y       yy     xx 


M    =  (l-v)  DW, 
xy  xy 


As  before,  the  displacements  are  expanded  in  Fourier  series 
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U(x,9)  =   Z   Un(x)  cos  n6  (A. 7a) 

n=0 


V(x,8)  =   Z   V^x)  sin  n9  (A. 7b) 

n=l 


W(x,0)  =   Z   W^x)  cos  n6  (A. 7c) 

n=0 


The  terms  comprising  the  integrand  of  (A. 4)  can  be 
expressed  in  terms  of  the  displacements  by  means  of  (A. 5) 
and  (A. 6),  leading  to 

2  W 

N  e   =  B[U,   +  v(V,   -  -)  U,  ] 

xx      x      y   a    x 

w  ^         w 

N  e   =  B[  (V,   -  -)   +  v(V,   -  -)  U,  ] 
y  y        'y    a         Y        a'   'xJ 

N   y      ^r^B  (U,   +  V,  )2 
xy '  xy     2      y     x 


My   =  -D  (W,2   +  vW,   W,   ) 
xAx         'xx       xx   yy 


M  Y   =  -D  (W2.    +  vW,   W,   ) 
yAy         yy      yy   xx 


M   x      (1-v)  DW2 
xyAxy  xy 


(A. 8) 


Replacing  the  displacements  in  (A. 8)  by  the  series  of 
(A. 7)  yields 


00      oo 


N  e   =  B  I    Z  [U,nU,m  +  ±    (n/  -  Wn)  U,m] 
n=0  m=0    X   X    a 


cos  n9   cos  m0  (A. 9a) 
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00  00 


N   e      =   B    Z         I    [i    (nV11   -   W11)  (mV™   -   W™)    + 

y  y        „_ri  m_n  =^ 


n=0   m=0    a' 


-  U,n    (mVm  -   W™) ]    cos   n6    cos   m6  (A. 9b) 

a        x 


oo  oo 

N      Y         =    ±Z±  b    Z  Z     (V,n    -    S.  u11)  (V,m    -   5L  Um)  • 

xy '  xy  2  t         i         x        a  x        a  , 

J      J  n=l   m=l 

sin   n0    sin  m0  (A. 9c) 

00  00 

M   x      =    "D    Z         Z     (W,n   W,m      -    4-  n2wV?      )• 
x   x  n=0   m=0        xx      xx        a2 

cos   n0    cos   m0  (A.9d) 

M    X      =    "D    Z  Z     (*^£)    wV1    -    v4  WnwT       )  ' 

y   y  n=0   m=0      a4  a2 

cos    n0    cos   m0  (A.9e) 

00  oo 

MX         =     (1-v)    D    Z  Z       — =■  W,    W, 

xyAxy  ,         ,       2       'x    'x 

J      J  n=l  m=l   a 

sin   n0    sin  m0  (A.9f) 


The  potential  energy  of  each  mode  is  obtained  by  substitut- 
ing (A. 7c)  and  (A. 9)  into  (A. 4)  and  performing  the  circum- 
ferential integration.   When  m,  n  ^  0,  the  potential  energy 

4-  V* 

of    the   n        mode    reduces    to 
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nn  =  ™  /L  (B[(U,n)2  4-  J  (nVn  -  Wn)  0  *]  + 


Z.  X      a.  _  *• 


*U.x  <»V»-W»>]  +D[(W,«x,2-^n2w7xxl 

2       2 
+  2(l-v)  D  ^  (W,£)  (A. 10) 

a 


4      2     2 
+  D[-j  (W  )  -  V-j  W  W,   ]  }  dx 
a  a 


When  m  =  n  =  0 ,  the  potential  energy  is 


n°  =  w  {bku ,0)2   -  ^w°u  o  +  J_  (w0)2] 

0        x     a     x    a2 
+  D(W,°  )   +  2PW°)   dx  (A. 11) 

.X.  .X. 


The  total  potential  energy  is  given  by 


n  =  z  nn 

n=0 


The  principle  of  Minimum  Potential  Energy  applies  to  each 
mode  individually  as  well  as  to  the  collective  sum  of  the 
potential  energies  of  all  the  modes.   Furthermore,  the 
potential  energy  of  a  mode  may  be  multiplied  by  an  arbitrary 
constant  without  changing  the  conditions  necessary  for  its 
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minimization.   Accordingly,  n   is  multiplied  by  one-half. 
This  allows  us  to  write  a  general  expression  for  the 
potential  energy  which  is  valid  for  all  values  of  n.   For 
convenience  the  terms  are  grouped  into  two  integrals, 
leading  to 


nn   =    ™    fL    {B[(U?)2    +    *£    (nVn    -    W11)    U    " 

C.  -.  X  CL  X 

+   J_    (nVn    -    W11)2]    +    D[(W,n    )2    -    2\  nVW/n 
2  xx  2  xx 

a  a 

4  2 

+   ~    (VP)     ]    +    26QnPWn    }    dx  (A. 12) 

a 

+   ^/L    {±^B<V«   -   ^Un)2 
z    0  I  x        a 

n2  n    2 

+    2(l-v)    D  2j    (W,£)     }      dx 

a 

where    6n      =1  when   n   =    0    and    6rt      =0   when   n   ¥■    0. 
On  On  7 

Let   II      be   defined   by 

ffn   =   -2—  nn  (A. 13) 

iraB 

Expressing  (A. 12)  in  terms  of  the  nondimensional 

variables  of  (2.21)  and  using  (A. 13)  give 

l  2  ? 

?rn         rr/n\  o/n  n  >         n         /      n  n>^ 

II      =  /      {  (\itfr )      +    2v  (nv      -   w    )    u,_    +    (nv      -   w    ) 

0  K  K 

2  0     -n         rt  A         n      ^ 


+   —  a2     [(w,£r)  2vn2wnw,^r    +   n4(wn)     ]+    26n^pwn}    d? 

(A.14) 


KK,  *v«    w    w,^    -r    xt     vw    ,     JT    ,U()nJ 


,      1-v     r      n  n.,2         1  2    2-      n    2 

+   /      ~2~    [v,?    -   nu    )  ]      +  £-    (1-v)    a   n    (w^)    }    d£ 
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Conventional  differences  will  now  be  used  in  conjunc- 
tion with  equation  (A. 14)  to  derive  a  set  of  difference 
equations.   In  order  to  discretise  the  potential  energy 
integral,  let  the  £  -  coordinate  of  the  n   mode  be  divided 
into  M  elements  each  of  length  d.   At  the  center  of  each 
element,  locate  a  node  point  K.   Then  II   may  be  expressed 
by  the  following  series. 

M  2 

rj-n        T.         „       ,      n .  «     ,      n  n  >  ,      n  x 

II      =   Lim      I       (u,r)      +   2v  (nv„   -   wv)(u,r)v. 

M->oo    K=l  K 

/      n  n.  ^  i       2r/      n      2         0i    2    n  ,      n    > 

+    (nvK   -   wK)      +na    [<w,u)K   -    2vn  wR(w,^)K 

+   n4(w£)2]    +    260npRw^   +   ^    [(v,^)K   -  (A. 15) 

n,2         1  ..       ,       2    2,      n    2 
nuR]       +    g-d-v)    a    n       (w'r)K-r    d 

where   d   =    i/M.      The   potential   energy   IT     can   be    approximated 
by   taking   M  to   be    finite    and   using   central    differences    to 
express    the   derivatives. 

The   principle   of   minimum  potential   energy   demands 

.,.,,-,.,  nnnnnn  nnn 

that    the   displacements    u,  ,    v,  ,    w,  ,    u„  ,    v„  ,    w ,,,...,  u.. ,    v..,    w.. 

1    1    1    2    2    2       M   M   M 

minimize  the  potential  energy  TT  ,  i.e. 


xtt^           v  311  ,    n  311       -    n  3U  ~    n        n                   ,,    ,  ^  N 

611      =      I      6uT,   +   6vT.   +   6wTr   =    0                   (A.  16) 

«._-,  *  n  K  .  n    K  -  n  K 

K=l  3uR  3vR  3wR 

Since  the  virtual  displacements  6u, ,6v1 , 6w, , . . . , 6u  , 6v  ,  6w^ 
are  arbitrary,  it  follows  that 
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^   =0 


3(UK) 


911    =  0  (A. 17) 


»s  /  n. 
9(vK) 


=0  K  =1,,.. ,M 


8(w£) 


When  the  derivatives  in  (A. 15)  are  replaced  by  their  con- 
ventional central  difference  approximations  (A. 2),  the 

potential  energy  becomes 

n      n  ~  n    n    ~ 

Fn    ,       ,UK+2  "  UKN2  ,UK   -    UK-2,2    ... 

n  ={•••+  ( — )  +  •••  +  ( — )  + 

2d  2d 

un    -  un 
.  ~  .   n      n   .  ,  K+2    K.    _  ,   n    n.  . 
+2v(nvK+1  -  wR+1)  ( — )  +  2v(nvK  -  wR) 

un    -  un                       un  -  un 
,  K+l     K-l,    o  /   n       n   WUK     K-2N 
( ^ )  +  2v(nvK_1  -  wK.1)( jg )  + 

,   n     n  2        i   2r/WK+2  "  2wK+l  +  WK,  2 
+  (nvK  -  wK)   +  •  •  •  _  a  [  ( -2 ) 

i    ,WK+1  -  2wK  +  WK-1,2  :   #wg-  2wg_,  +  wg_   2 

+  ( 5 )   +  ( 5 )  + 

d  dz 

,       2      n          ,WK+2    '    2WK+1    +   WKV     ,       2    n     ,wK+r2wK  +  WK-1N 
-2vn      wR+1     ( - )-2vn  wR    ( ) 

d  d 

-.  w£   -    2w£    ,    +  w£    ~  .  ~ 

~      2    n  ,    K  K-l  K-2.,  4,    n»2,,         ,  or  n 

"2vn  Vl    ( ~2 >  +  v.;-+n    (w    )     l  +  ---+2%nPKw 

d 

vn    „  -  vn  o  vn    ,  -  vn -  -  o 

J.1-V     r  /    K+2         K  n         2  /K+l         K-l  nN2 

+  ,**+T-    [( 2d nUK+l}       +     ( 2d nUK}       + 

n        n  n  n 

,VK~VK-2             n         2                i                2    2     r  ,WK+2~  WK,2^         . 
( — nuK_1)      +...+Ffl-v)a  n      [( )+...+ 

n  n 

(    *    0,    *    ^)     ]    +    ...}    d  (A. 18) 
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Only  the  terms  containing  u„,  vv,    and  w   have  been  displayed 
since  the  displacements  at  the  other  stations,  K+l,  K-l, 
etc.,  do  not  enter  into  the  variation  of  u.,,  vv,    and  w„. 
The  axial,  circumferential,  and  radial  difference  equations 

are  obtained  by   equating   to   zero   -Kj—^r    ^rrry    ^f^*      Tlie   result- 

(uk)   (vi3  1wk) 
ing  equations  are 

Axial: 

n  n  n  n  n 

UK+2    ~    2UK    +    UK-2  14-v       ,VK+1    "    Vlt 

5 +  -~-  n( )        (A. 19a) 

4d  2d 

n      n 

1-v   2  n     ,WK+1  "  WK-1,  n 
— ^—  n  uv    -  v( )=0 

2      K  2d 

Circumferential : 

n      n  n      _  n    n 

2  n    1  +  v   ,Vl_Vlu  1-v  ,VK+2  "  2vK    VK-2X 
-n  v   -  -*—  n( )+  — —  ( - ) 

K     2         2d  2  4d2 


+  n2wn  =  0  (A.  19  b) 


Radial: 


n      „  n      ,  n    ,  n       n 
1_  a2    WK+2  -  4WK+1  +  6WK  -  4WK+1  +  WK-2. 
12  d4  ) 

n      ~  n     n 

,  2  ,WK+2  "  2WK  +  wK-2„.   4   n 

-2n   ( )+  n  wv  + 

4d2  K 

n      .n      ,  n    „  n      n 

v   2^2  ,WK+2  "  4WK+1  +  6WK  "  4WK-1  +  V2a_ 
-  n  d   ( )]- 

z  dq 

n      n 

x         n     n     n    ,,.UK+1    UK-1, 

o„  p   +  nv   -  w   -  w  +  v( ) 

On^K      K     K     K     v     2d 


(A. 19c) 
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Let  us  compare  the  above  difference  equations  with 
(A. 3)  which  were  derived  directly  from  the  differential 
equations.   Note  that  the  second  derivatives  in  equations 
(A. 19)  are  represented  by  expressions  of  the  form 

(f'£5}K  =  772  (fK+2  "  2fK  +  fK-2) 
4d 

as  opposed  to  the  conventional  central  differences  of  the 
second  derivatives  in  (A. 3). 

In  addition,  (A. 19c)  contains  a  term  which  has  no 
counterpart  in  (A. 3c)  namely 

v   2^2  2  ,wK+2  "  4WK+1  +  6WK  "  4wK-l  +  WK-2N 
__  a  d  n   ( ) 

2 
Were  it  not  for  the  pressence  of  n  ,  this  term  might  be 

discarded  as  a  higher  order  term.   However,  the  problem 
under  investigation  encompasses  response  modes  having 
several  circumferential  waves;  and  the  term  may  not  be 
small  compared  to  other  quantities  which  appear  in  the 
radial  equilibrium  equation. 

Consider  now  the  modified  central  differences.   A 
set  of  difference  equations  will  be  derived  by  minimizing 
the  potential  energy  in  modified  difference  form.   Conse- 
quently, a  series  representation  of  (A. 14)  must  replace 
the  integral  formulation.   The  right  side  of  (A. 14)  was 
written  as  the  sum  of  two  integrals  for  reasons  which  will 
now  become  apparent.   In  order  to  obtain  a  series  which  is 


194 


equivalent  to  the  first  integral,  let  the  ^-coordinate  be 
subdivided  such  that  a  whole  station  node  is  located  at 
the  center  of  each  segment.   This  procedure  is  schematic- 
ally illustrated  in  Figure  A. la.   In  approximating  the 
first  integral,  the  half  station  at  each  end  is  neglected. 
Doing  so  does  not  affect  the  difference  equations  that  must 
hold  at  the  interior  points.   Hence,  the  first  integral 
can  be  written  in  the  form 

M-l        2  2 

1st  Integral  =   Z   {(u'r)K  +  2v (nvR  -  wK)(u'r)R  +  (nVK~WK^ 

K=  2. 

1   2,.   n  .2       2   n  ,   n  .      4.  n.^    ox      n,  , 
+  12  a  [(W'«JK  "  2Vn   WK  (W'^>K  +  n  (WK}  ]+  260nPKWK}  d 

The  second  integral  of  (A. 14)  is  approximated  by  using 
finite  elements  which  are  centered  at  half  station  nodes, 
as  shown  in  Figure  A. lb.   The  series  approximation  of  this 
integral  is  given  by 

o  j  -r       t     r  /  1-v  ,  ,   n.z      n,z 

2nd  Integral  =   I  I— —  [(v,-)   -  nu  ] 

K=l/2,3/2,. ..       l  ? 

1        2  2    n  2 
+  -  (1-v)  a  n  (w,t)KJ  d 

Combining  these  two  series  leads  to  the  equation  for  the 
potential  energy.   Replacing  the  derivatives  by  modified 
central  differences  of  (3.16)  and  minimizing  the  potential 
energy  lead  to  the  governing  difference  equations 
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4  M-3 


M-2 

— o — 


M-l 


M- 


TTTTTTT 


ELEMENT   #2 


FIGURE    A. la 
"WHOLE"     STATION     MESH- MODIFIED    DIFFERENCES 


1/2 


& 


.*1 


5/2 


ELEMENT    4*=  2 


M-l/2 

M 


7h 


FIGURE    A. lb 
"HALF"  STATION    MESH- MODIFIED     DIFFERENCES 
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„  n        n     .    1-v    1  ,  n      n. 
-I<UK+3/2"    K+1/2  "  UK+1/2     ~  "  d  <VK+1  "  V 

1-v      2  n  1/H        .  ik       . 

— T  n    UK+1/2"V    d(WK+rWK)=    ° 


2  1      n        1-v        1,    n  n  x,  1-v     1„.    n  ~    n  n      4 

"n    d   VK    +    T"  n    d(UK+l/2-    UK-l/2)+—    d2(VK+l"   2VK    +    VK-1> 

+   nw"   =    0  (A. 20) 


1      2  ,  1     ,   n  .n  ,    n   _    „    n  n      » 

12    a    [74     (WK+2    "    4WK+1    +    6WK         4WK-1    +    WK-2> 
d 


o    2     1      /    n  »n,riv  4    n.  , 

-2n      3    (wR+1    -    2wR   +   w^)    +   n   wR)  ] 

d 


=  60nPK  "  nVK  "  WK  +V  |  K+l/2  "  UK-l/2) 


These  equations  are  identical  to  the  static  version 
of  (3.17),  equations  obtained  by  applying  the  modified 
differences  to  the  governing  differential  equations  (3.15) 

When  the  modified  differences  are  used  to  approxi- 
mate the  derivatives,  the  set  of  difference  equations 
obtained  by  satisfying  the  differential  equations  at  a 
finite  number  of  points  are  the  same  as  the  difference 
equations  obtained  by  satisfying  the  condition  of  minimum 
potential  energy.   On  the  other  hand,  when  conventional 
central  differences  are  used,  the  difference  equations 
obtained  from  the  differential  equations  do  not  agree  with 
those  obtained  by  the  energy  method.   Consequently,  the 
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modified  differences  should  yield  more  accurate  solutions 
than  conventional  differences  for  the  same  mesh  length. 
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APPENDIX  B 

ADAPTATION  OF  MODIFIED  DIFFERENCES 
TO  THE  NONLINEAR  EQUATIONS 

I.   EVALUATION  OF  NONLINEAR  TERMS 


In-plane  Forces,  \i,  and  A2 


The  in-plane  forces  were  expressed  in  Chapter  III 
by  (3.5).   Let  n_  and  n   be  defined  as  whole  station 
quantities,  and  let  nr„  be  a  half  station  quantity. 
Replacing  the  derivatives  in  (3.5)  with  the  appropriate 
modified  difference  expressions  from  (3.18)  leads  to 

,  n.     1  /  n    n   .     ,   n    n. 
(Vk  =  d(Uk  '  Uk-1>  +  V(nVK  "  WK} 


(B.la) 


+  4  +  CITT   +  *CC*   +  C;TT  ) 

J\  K.  K  K 


.    n.       n    n    v  ,  n    n   v 
(Vk  =    nVK  "  WK  +  d  (Uk  "  Uk-1} 


J*         +  Cn     +  v(Cn    +  Cn    ) 
'TT__     ITT^       XXV  IXX_/ 

K  K.  K.  i\ 


,    n  .     1-v  rl  ,  n       n.      n    ^n 

(lVk  =  —   [d  (vk+i  "  vk)  -  nuk  +  cxTk 

+  cn  +  cn   ] 

k        k 


(B.lb) 


(B.lc) 


The  nonlinear  terms  in  the  axial  and  circumferential 
equations  are  contained  in   XI   and  X2    ,    defined  by  (3.8) 
and  (3.9).   Equation  (3.20)  requires  that  Aln  be  defined 
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at  half  stations,  and  equation  (3.21)  requires  that 
X2      be  defined  at  whple  stations.   Approximating  (3.8) 
at  the  half  station  k  and  (3.9)  at  the  whole  station  K 
results  in 


\-\n   =  L(rn  -    rn      \    +  —(rn  rn         \    + 

xk    dv  XX^,    XX./  d^IXXVJ,    IXXT/ 

K.+  X       J\  K+X         K 

v[i(C^     -  CL,  )  +  ^-(C1}      -c2m  )]      (B.2a) 
K+l      K  K+l     K 

1-v  .  _n      _n      _,n    .   v  ,  n      nx 

+  -2-n(CXTv  +  CIXT  +  CITXJ-  d(WK+l"  V 

k  k       k 


A2K=  "n[C?T,  +  CITT+U(CSx/CIXX  »1  + 

J\  i\        K       K 

T1  'd  (C5Tk-C^k-l)+  d  ^XT^IXT^'   <B-2b» 

+  3(CITX  "  CITX,   »  1  +  ™K 
k      k-1 


Examination  of  (B.l)  and  (B.2)  shows  that  C™ ,  CTT_, 
Cvv,  and  C    must  be  defined  at  whole  stations  while  C    , 

XX         XXX  XXI 

Cvm,  and  CTm   must  be  defined  at  half  stations.   As  seen 

XI  XIX 

from  (3.4a)  and  (3.4d),  C     and  C     are  obtained  by 

XX  XXX 

summing  appropriate  products  of  the  form 

i    j 
w,  w  ,J 

and 


—  i    D 
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Since  C^v  and  c"    are  defined  at  a  whole  station,  w,r 

XX        XXX  ^ 

and  w7  must  be  calculated  at  whole  stations.   However,  w 
is  defined  at  whole  stations  only;  and  the  modified 
difference  form  defines  w,   at  half  stations,  not  whole 
stations.   A  similar  situation  exists  for  the  half  station 
nonlinear  terms  C    and  CIXT,  terms  which  can  be  evaluated 
only  if  w   is  known  at  half  stations.   This  creates  some 
difficulty  when  utilizing  the  modified  finite  differences 
to  solve  the  nonlinear  equations  numerically.   The  diffi- 
culty is  overcome  through  the  use  of  interpolation  formulas 
or  simple  averaging. 

Consider  the  problem  of  determining  w,r  at  a  whole 
station.   With  reference  to  Figure  B.l,  let  w,  _,  and  w, 
be  the  radial  displacements  at  half  stations  k-1  and  k, 
respectively.   The  required  derivative  of  w  ,  which  we 
shall  denote  (w,  )  ,  is  given  by 

(wVk  ■  d(wk  -  wk-i>  <B-3) 

By  expanding  w   in  Taylor  series,  the  quantity  (w£  -  w£_, ) 
can  be  expressed  as  a  linear  combination  of  w   0,  w„  ,  , 
w   ,  ,  and  w   _.   Using  the  Taylor  series  to  express  w„,  we 
have 

2 

n    n   d,  n  .     1  ,d.  ,    n   v     „  ,,3>        ,_,  .v 
wK  =  wk  -  j(w,e)k  +  j{j)     (w,^)k  +  0(d  )        (B.4) 
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FIGURE  B.l 
FINITE  DIFFERENCE  STATIONS 


where  (w,,,).  and  (w ,,.,.).  are  the  first  and  second  deriva- 
€  k         ££  k 

tives  of  w   at  half  station  k.   The  notation  0 (d  )  means 

3 
that  the  truncated  terms  are  of  order  d  .   The  second 

derivative  (w,_p),  can  be  expressed  in  terms  of  (w,--)   or 
(W^?}K-H1  aS 

(w*Wk  =  (W^C}K  +  f  (W^C^^}K  +  0(d2)  (B'5a) 

(W^^k  =  (W^>K+1  "  |^5«}K  +  °(d2)  (B-5b) 


Adding  the  above  two  equations  yields 


(w"EC)k  =  |[(«"SC)K+  (W"«'K+1]  +  0(d2) 


(wVJ^  -  w^,  -  wv  +  w^  , )  +  0 (d  )  (B.6) 


^2   K+2     K+l     K     K-l 
2d 


Also,  we  have  that 


<WVv     =    5<WK+1    "    "k»    +    °<d2)  <B-7) 
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Substituting  (B.6)  and  (B.7)  into  (B.4)  leads  to 

n    9  .    n    n   ,    1  .    n      n      .,,3.    .-,   0. 
Wk  =  I6(WK  +  WK+1>  "  I6(WK-1  +  WK+2)  +  ° (d  >    (B'8) 


It  follows  immediately  from  (B.8)  that 


n      9  ,  n      nv    1  ,    n      n   .    n  ,  ,3,  ,_,  n . 
Wk-1  =  16(WK-1  +  WK}  "  T6(WK+1  +  WK-2>  +  ° (d  )  (B'9) 


Finally,  the  interpolation  formula  for  the  derivative  is 
obtained  by  substituting  (B.8)  and  (B.9)  into  (B.3), 
resulting  in 


(WVk   ■    ira[10(WK+rWK»     ".    K+2-^-2'1    +    0(d2»     <B-10a) 

When  it  is  necessary  to  calculate  w   at  a  half  station, 

either  (B.8)  or  (B.9)  may  be  used. 

At  the  boundaries  the  first  derivative  of  w  may  be 

approximated  by  using  a  forward  or  backward  difference  with 

(B.8)  or  (B.9).   Since  w,  =  w„  =  0,  the  result  is 

l  M 

(wVk=1   =    <!>  [I6  W2   *   T6iWo  +  "31    +   0(d)  <B-10b) 


and 


<W"C>K=M   =    "     (l>  [R  WM-1    "    T?(WM+1    +   WM-2>+    °<d)  (B-10C) 


The  procedures  for  evaluating  the  Fourier  coefficients 
of  the  nonlinear  terms  defined  by  equations  (3.4a)  -  (3.4g) 
are  described  in  Appendix  B  of  Reference  [30].   For  example, 
Cxx  at  whole  station  K  is  given  by 
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/^.n      l      n  x  r  c,  i+n.      c,    i-n  . 

(CXX)K  =  4  Jo^'^K^  (W'C  )fk  +  y  (W4     }K]   (B.ll) 


where 


c 

n  = 


1  if  n  >  0 

0  if  n  =  0 

1  if  i  f   n 

2  if  i  =.  n 


The  terms  of  the  series  can  be  evaluated  by  utilizing 
(3.18),  (B.8),  (B.9),  or  (B.10).   The  other  nonlinear  terms 

CIXX  '  C?T  '  C?TT  '  Cxt  '  C?XT  '  and  C~ITXV  Can  be  evaluated 
from  equations  similar  to  (B.ll). 

n    n         n 
Nonlinear  Terms n£  ,  ^£9  ,  and  n6 

The  nonlinear  terms  in  the  radial  equation,  (3.22), 

are  embodied  in  the  terms  nQ  ,  n-  ,  r\r~    ,  and  r\a    .   Equation 

9K    CK    C9K       6K 
(B.2b)  is  used  to  compute  the  first  of  these  terms.   The  non- 
linear terms  ric  and  iv  are  defined  by  (3.4h)  and  (3.4j). 

The  calculation  of  (n,.)   involves  the  summation  of  products 

^  K 
of  the  form 


<n£)R  I(w3u)K  +  (w?u)K] 


The  coefficients  C^xx  ,  C^TT  ,  C^XT,  and  C^TX  are 

not  truly  nonlinear  since  they  are  sums  of  products  that 
contain  derivatives  of  the  initial  imperfection.   However 
they  are  computed  in  the  same  manner  as  the  nonlinear 
terms . 
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for  appropriate  values  of  i  and  j,  and  the  calculation  of 

(ru)   requires  summing  products  of  the  form 
9  K 

(ni)   [-j2(wj  +  wJ)  ] 
y  K  K 

The  summations  are  defined  by  expressions  similar  to 

(B.ll).   The  terms  in  the  products  can  be  computed  from 

(B.la),  (B.lc),  and  the  difference  approximations  of 

(3.18)  . 

The  remaining  nonlinear  term  (nro)   requires  the  use 

t,v   K 

of  a  simple  average.   Equation  (3.4i)   expresses  (nrQ)   as 

£6  K 

oo  OO  00 

(  I    2n^Q   sin  ie)[  Z  -  j  (w^  +w?  )  sinj  6]-  Z  n"    cos  ne 
i=l   C9K         j  =  l      ?K  k  n=0  ^6K 

(3.4i) 
However,  (B.lc)  defines  n™  as  a  half  station  quantity.  We 
take  rira   as  tne  average  of  nrQ   and  nro    ;  thus,  (3.4i) 

ceK  cek     eek_x 

can  be  written  as 


oo  oo  oo 

Z  rirn  cos  n6  =  [  Z  (n^Q  +  n^0    )sin  i6][  Z  -j(w?c  H-w^   ) 

n    K6K  i=l   ?9k    ?9k-l  j  =  l      *>K    ?K 

n=0  J 

sin  j6]     (B.12) 

where,  as  before,  w,    is  the  first  derivative  of  w   at 
whole  station  K  and  can  be  evaluated  by  using  equation 
(B.10).   From  (B.lc)  it  follows  that 

n      n       1-v  r  n      n       .  n    n   \       lo    -i^\ 

n^ek  +  ncek_1  =  1"  [vk+i  "  vk-i  "  n(uk  +  uk-i>  (B'13) 

CXTk    CXTk_1+  CIXTk+  ^XT^   CITXk  CITXk.1J 
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Finally  (rirfl)   can  be  evaluated  by  using  a  formula  from 
t«  K 

Appendix  B  of  Reference  [30],  similar  to  (B.ll). 

II.   Boundary  Conditions  on  Displacements 

Consider  the  case  in  which  the  axial  displacements 
are  zero  at  the  ends.   Since  un  is  not  defined  at  the  whole 
stations  K=l  and  K=M  at  the  ends  of  the  shell,  interpolation 
formulae  are  required.   The  finite  difference  stations  near 
the  end  £  =  0  are  shown  in  Figure  B.2.   The  displacement 
u,  _.  can  be  expressed  by  the  Taylor  series 

n  n  d  ,   n    >  1  ,d»     ,    n      x  .  ,,2,  /t-.    -i  <  \ 

Uk=l=    UK=1+    2(U'C}K=1+    2(2}     (U'^}K=1+    °(d    >  (B*14) 


Now 


n  n 

u?c  =      k=1    „ lizi+0(d2)  (B.15) 

?K=1         d 


and  using  forward  differences  gives 

(u")     =  C  K"2 Ozl+0(d)  (B.16) 

^    K=l  d 


We  also  have 

n      n 

n        uk=2  "  uk=l      2 
(un      =  JLJ ^L+0(d2)  (B.17) 

*  K=2        d 


To  obtain  the  desired  boundary  condition,  (B.15)  and  (B.17) 
are  substituted  into  (B.16),  (B.15)  and  (B.16)  are  sub- 
stituted into  (B.14),  and  (u  )K=1  is  equated  to  zero.   The 
result  is 
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K=l  K=2 


FIGURE  B.2 

FINITE  DIFFERENCE  STATIONS 
AT  THE  END  £  =  0 


Uo  +  2U1  "  3  U2  +  0(d3)  =  °  (B.18) 


Similarly,  the  condition  that  u  =  0  at  £  =  i  can  be  ex- 


pressed by 


u"  +  2i£   -  i  u"    +  0(d3)  =  0  (B.19) 

m     m— 1    3   iti-z 


When  the  shell  is  clamped  at  K=l  and  K=M,  the  slopes, 
(w'r)^_i  an^  (w i  r)  is-™   roust  be  zero.   Equating  (B.lOb)  and 


(B.lOc)  to  zero  yields 


wn  -  9w"  +  w*  =  0  (B.20) 

O      2     3 


WM  0     ~     9WM  1  +  WM  1  =  °  (B.21) 

M-2      M-l     M-l 
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APPENDIX  C 


DEFINITIONS  OF  CONSTANTS  AND  MATRICES 


In  equation  (3.22) 


Nr 
dl  "  X  +   *   6KH.C2 


d9  =  -  £ 

2     a 


d3  "  "  12  d4 


.3      1  a2  ,  2    2  , 


2 

d.  =  -  ~  2L.  (6  +  4n2d2  +  n4d4) 
5      12  d4 


In  equation  (3.23) 


2 

dfi  =  -  ~  [C.e  (2-en"1)  -  C,] 
6      d2    4  5 


d_  =  -  C„n2  -  C,  (n2-l)2  -  C, (en2-  l)2  -  2d, 
7       4       1  3  6 

d8  "  "  TT  [C4(1-en2>  +  C5! 


2 

dg  =  nC3  (1-en  ) 


In  equation  (4.2) 
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An  = 


dn(— ) 


1-v 


B 


Nr 

-2  -  d2n2(^)  +  I    6RH  d2n2(C4n2+C5) 

j=l    j 


,1  +  V\   ^ 

(— )  dn 


Nr 

-l(l-v)  -  d2n2  +  I    6V„   d2n2C, 

j=l  KHj      3 


Cn  = 


dn(— 3 


l-v 


,n 


In  equation  (4.8),  P,  is  given  by  the  following  expressions 


Simply  supported  ends: 


*1  = 


-1 


,    ,2  2  /l~v v 
1  +  d  n  ( — -) 
2 


j  / 1+  v  \ 
-dn  (-2—) 


,    ,2  2 ,1-v. 
1  +  d  n  (-r-) 


Clamped  ends: 


*?- 


-4/3 


.    ,2  2  ,1-v, 
4  +  d  n  (— — ) 


,    ,2    2,1-Vv 
4  +  d  n  (-y-) 
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In  equation  (4.11)  matrices  H   and  0   are  given  by 
the  following  expressions. 


Simply  supported  ends : 


Hn  = 


.    ,2  2  ,1-v. 
1  +  d  n  (-2") 


dn(— ) 


dn(-y-) 


2  2 
1-v  +  d  n 


0n  = 


dn(-2~ ) 


1-v 


n 
rM-l 


f  =  g,.  , 


d<cxx  +  c?xx> 


M 


Clamped  ends : 


Hn  = 


.    .2  2, 1-v. 
4  +  d  n  (-j-) 


dn(-j-) 


dn(~2— ) 


^.2  2 

1-v  +  d  n 
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on  = 


4/3 


dn(~2— ) 


1-v 


"K 


n 


In  equation  (4.12) 


Rn  = 


dn(-y-} 

l-v  T~l2~2 

— =—  +  d  n 


{    2    ' 


,1-v,  ,  -,2  2 

( )  +  d  n 

2 


n 
s   = 


A2 


n 
"m 


l-v    ,2  2 

— ^—  +  d  n 
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APPENDIX  D 

COMPUTER  PROGRAM 

I.   FORTRAN  VARIABLES 

Table  XI  lists  the  major  Fortran  variables  and  their 
algebraic  counterparts  from  the  main  body  of  the  disser- 
tation, if  such  a  counterpart  exists.   The  equation  in 
which  the  algebraic  counterpart  was  first  defined  is  indi- 
cated in  the  remarks  column.   The  more  important  program 
control  variables  are  also  briefly  described. 

II.   DESCRIPTION  OF  THE  PROGRAM 

The  computer  program  has  been  written  in  FORTRAN  IV 
language  for  operation  with  the  IBM  360  MODEL  67  Computer. 
The  program  was  compiled  by  the  FORTRAN  IV  G  Level  1, 
Mod  3  compiler.   The  program  is  dimensioned  to  accomodate 
twenty  equally  spaced  axial  stations  and  fifteen  arbitrary 
Fourier  modes.   The  program  will  handle  up  to  three 
indentical  ring  stiffeners,  located  so  as  to  divide  the 
shell  into  segments  of  equal  length.   The  memory  require- 
ments and  the  execution  time  are  contained  in  part  IV  of 
this  appendix. 

The  program  solves  a  set  of  three  nondimensional , 
nonlinear  equations  by  an  iterative  procedure.   The  depen- 
dent variables  are  expanded  in  Fourier  sine  or  cosine 
series  in  the  circumferential  coordinate.   The  nonlinear 
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terms  are  thus  converted  to  products  of  Fourier  series. 
The  products  are  expanded  to  obtain  a  single  series  for 
each  such  product.   Derivatives  with  respect  to  the  axial 
coordinate  are  evaluated  numerically  by  using  a  set  of 
modified  staggered  node,  central  differences. 

The  radial  equation  of  motion  for  each  mode  is  inte- 
grated in  the  time  domain  by  using  Newmark * s  beta-method. 
After  determining  the  radial  displacements,  the  nonlinear 
terms  in  the  axial  and  circumferential  equations  are 
evaluated.   The  axial  and  circumferential  displacements 
are  then  computed  by  Gauss  elimination.   Finally,  the 
radial  accelerations  are  computed  from  the  radial  equilib- 
rium equations;  and  a  convergence  test  is  made.   The 
process  is  repeated  until  convergence  is  achieved  or  until 
a  specified  maximum  number  of  iterations  have  been  com- 
pleted.  If  the  solution  fails  to  converge  in  the  maximum 
allowable  number  of  iterations,  the  time  increment  is 
multiplied  by  two-thirds  and  the  procedure  is  repeated. 
Whenever  convergence  occurs  in  less  than  some  specified 
minimum  number  of  iterations,  the  time  increment  is  in- 
creased by  15  per  cent. 

The  output  consists  of  the  maximum  circumferential 
stress  at  G  =  0   attained  during  the  solution  together 
with  the  time  and  station  at  which  the  maximum  stress 
occurred.   In  addition,  the  displacements  and  in-plane 
stress  resultants  are  printed  at  specified  intervals. 
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Main  Program 

The  main  program  performs  the  following  functions: 

1.  Controls  the  logical  connections  between  the 
subroutines . 

2.  Adjusts  the  length  of  the  time  increment. 

3.  Sums  the  Fourier  coefficients  in  the  series  for 
the  in-plane  stress  resultants. 

4.  Computes  the  circumferential  stress  due  to 
bending. 

5.  Computes  the  maximum  circumferential  stress 
together  with  the  location  and  time  of  its  occurrence. 

6.  Prints  out  the  solution  at  specified  time 
steps  and  after  completion  of  the  solution. 

Subroutine  START 

This  subroutine  performs  the  following  functions: 

1.  Reads  in  the  data  cards. 

2 .  Prints  the  problem  description  and  program  con- 
trol parameters. 

3.  Sets  the  initial  values  of  all  array  elements 
to  zero. 

4.  Calls  subroutines  WIDER  1  and  WIDER  2,  or 
WI1SIN  and  WI2SIN,  which  compute  and  store  derivatives  of 
the  initial  imperfection. 

5.  Sets  certain  program  control  parameters  to  proper 
initial  values. 
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Subroutine  WDEF 

This  subroutine  performs  the  time  integration  by 
Newmark ' s  beta-method. 

Subroutine  MODES 

Subroutine  MODES  was  taken  directly  from  Reference 
[30]  with  the  permission  of  its  author. 

In  MODES,  the  arrays  that  define  those  sets  of  indices 
that  combine  to  equal  each  value  of  n  in  the  problem  are 
determined.   They  are  IS,  JS ,  ID,  JD ,  I JS ,  MAXS ,  MAXD,  and 
MAXSY.   MODES  is  called  during  the  first  iteration  of  the 
first  time  step  and  during  each  subsequent  iteration  until 
the  number  of  Fourier  terms  reaches  a  specified  number. 
The  arrays  computed  here  are  used  as  control  parameters 
by  subroutines  LAM1N2  and  TNETA  in  evaluating  the  nonlinear 
terms . 

Subroutine  LAM1N2 

The  quantities  LAM1  and  LAM2 ,  the  terms  which  com- 
prise the  load  vector  in  the  Gauss  elimination,  are  computed 
in  Subroutine  LAM1N2 .   The  derivatives  of  the  radial 
displacement  that  appear  in  LAM1  and  LAM2  are  first 
evaluated.   A  multiplying  and  summing  procedure  is  then 
carried  out  to  evaluate  CXX,  CIXX,  CTT,  CITT,  CXT,  CIXT, 
and  CITX,  the  Fourier  coefficients  of  the  expanded  products 
of  series.   Arrays  IS,  JS,  ID,  JD,  US,  MAXS,  MAXD,  and 
MAXSY,  which  were  generated  in  MODES,  are  used  as  control 
parameters  in  the  multiplications  and  summations. 
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Subroutine  PMATRX 

PMATRX  serves  as  a  control  program  for  the  computation 
and  storage  of  the  matrices  P,  DEE,  and  DST,  matrices  used 
in  the  Gauss  elimination  solution  for  the  displacements 
u,  and  Vj,.   The  subroutine  is  called  by  the  main  program 
during  the  first  iteration  and  whenever  the  expansion  of 
the  products  of  series  introduces  new  Fourier  modes.   The 
boundary  conditions  at  the  first  axial  station  are  in- 
corporated into  the  submatrices  of  P  associated  with  the 
first  station.   Subroutines  ABC  and  PANDD  are  then  called 
to  compute  the  rest  of  the  P  matrix  and  the  matrices  DEE 
and  DST. 

Subroutine  ABC 

This  subroutine  computes  the  elements  of  the  A,  B, 
and  C  matrices  for  a  given  Fourier  number.   These  matrices 
are  used  in  subroutine  PANDD  to  compute  the  elements  of 
P,  DEE,  and  DST.   At  stations  having  rings,  the  matrix  BR 
is  used  in  place  of  B. 

Subroutine  PANDD 

The  matrices  P,  DEE,  and  DST  are  computed  by  PANDD. 
Subroutine  MATMAT  and  INVMAT  are  called  to  perform  ele- 
mentary matrix  operations. 

Subroutine  UVDEF 

Subroutine  UVDEF  computes  the  displacements  u,  and 
vv,  which  correspond  to  the  new  radial  displacement,  wl . 
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UVDEF  first  computes  the  new  vector  GEE,  then  performs 
row  operations  on  the  vector  GEE  to  generate  the  X  vectors, 
and  finally  computes  u,  and  vv   by  back  substitution.   Sub- 
routines  MATMAT,  INVMAT ,  and  MATVEC  are  called  for  ele- 
mentary matrix  operations. 


Subroutine  TNETA 


This  subroutine  performs  the  following  functions. 

1.  Computes  TX,  TXY,  and  TY ,  the  Fourier  coefficients 
of  the  in-plane  stress  resultants. 

2.  Computes  WWX,  WWXTH,  and  WWTH,  the  Fourier  co- 
efficients of  the  ^-curvature,  the  twist,  and  the  circum- 
ferential curvature,  respectively. 

3.  Carries  out  a  multiplication  and  summation  pro- 
cedure to  compute  ETAX,  ETATH ,  and  ETAXTH,  the  Fourier 
coefficients  of  the  nonlinear  terms  in  the  radial  equilib- 
rium equations.   The  arrays  IS,  JS,  ID,  JD,  I JS ,  MAXS , 
MAXD,  and  MAXSY ,  which  were  prepared  in  subroutine  MODES, 
are  used  as  control  parameters  in  the  multiplication  and 
summation  process. 

Subroutine  ACELR8 

This  subroutine  solves  the  radial  equilibrium  equation 
of  each  mode  for  the  radial  acceleration  WDD0T1.   The  con- 
vergence test  is  performed  as  the  new  accelerations  are 
computed.   The  parameter  NCONV  is  set  to  zero  when  the  solu- 
tion   converges,  or  to  one,  if  the  convergence  test  fails. 
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The  program  stops  checking  for  convergence  as  soon  as  a 
test  fails. 

Subroutines  INVMAT ,  MATMAT ,  and  MATVEC 

These  subroutines  perform  elementary  matrix  opera- 
tions on  2  x  2  matrices  and  2x1  vectors.   INVMAT  inverts 
a  matrix,  MATMAT  forms  the  product  of  two  matrices,  and 
MATVEC  multiplies  a  vector  by  a  matrix. 

Subroutines  WIDER  1  and  WIDER  2 

These  subroutines  compute  and  store  the  derivatives 
of  the  initial  imperfection  required  in  subroutines 
LAM1N2  and  TNETA.   The  subroutines  are  based  on  the  assump- 
tions that  the  imperfection  in  mode  n  is  equal  to 

win  {1  -  cos  [2tt(N   +  l)g/l  ]}  cos  n8 

and  that  the  rings  are  located  such  that  the  shell  is 
divided  into  N   +1  segments  of  equal  length.   This 
means  that  the  imperfection  and  its  first  derivatives 
are  zero  at  the  ring  stations  and  do  not  produce  any 
strains  in  the  rings,  as  assumed  in  the  analysis.   Accord- 
ingly, there  is  a  finite  number  of  compatible  ring  sta- 
tions.  They  are  tabulated  in  the  program  listing  of 
subroutine  WIDER1.   The  program  listed  calls  these 
routines.   If  another  form  of  the  imperfection  is  desired, 
a  separate  routine  must  be  written. 
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Subroutines  WI1SIN  and  WI2SIN 

These  subroutines  are  substitutes  for  WIDER1  and 
WIDER2  and  were  used  in  computing  the  solutions  pre- 
sented in  Chapter  V.   These  subroutines  are  based  on  the 
assumption  that  the  initial  imperfection  in  mode  n  is 
given  by 

wi   sin(7T^/i)  cos  n0 

Subroutine  PRLOAD 

This  subroutine  computes  PLOAD  at  each  station  for 
a  uniform  time  decaying  external  pressure,  i.e. 

PLOAD (K)  =  PO  exp  (-RT) 

The  subroutine  may  be  replaced  by  the  user  for  other 
pressure  time  histories.   The  variables  PO,  R,  S,  and 
OMEGA,  in  labeled  common  block  Al ,  are  read  from  data 
card  number  ten.   These,  variables  are  not  used  anywhere 
else  in  the  program  and  are  available  to  the  user  in 
developing  an  alternate  subroutine  for  the  loading. 

Subroutine  AMPMAX 

Subroutine  AMPMAX  is  not  an  integral  part  of  the 
program.   It  was,  however,  employed  in  computing  the 
solutions  presented  in  Chapter  V.   It  should  be  used  only 
in  conjunction  with  WI1SIN  and  WI2SIN.   AMPMAX  computes 
the  amplification  of  those  Fourier  modes  for  which  an 
initial  imperfection  was  specified.   If  the  magnitude  of 
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the  initial  imperfection  for  some  mode  is  set  identically 
to  zero  by  the  input  data,  the  subroutine  will  compute 
an  amplification  of  zero  even  though  the  response  of  that 
mode  may  not  be  zero.   The  amplification  is  computed  at 
the  midspan  of  the  shell  or  at  the  whole  station  nearest 
to  the  modspan.   In  computing  the  amplification/  the 
initial  imperfection  is  assumed  to  be  equal  to 

wi   sm(— -)  cos  n0 

III.   INPUT  REQUIREMENTS 

The  layout  of  each  data  card  is  shown  in  Figure 
D.l.   The  first  three  cards  will  be  referred  to  as 
TITLE  1,  2,  3;  and  the  ramaining  cards,  as  DATA  1  -  13. 

The  program  uses  nondimensional  variables  throughout 
and  the  user  must  supply  the  data  in  nondimensional  form. 
The  contents  of  each  card  are  listed  below. 

Data  Cards 

TITLE  1-3:      The  description  of  the  problem  is  placed 

on  these  three  cards  and  will  be  printed 

at  top  of  the  first  page  of  output. 
DATA  1: 

MM  -  The  number  of  axial  stations ,  including 

the  stations  at  the  ends  of  the  shell; 

MM  <  20 
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NR 
INMAX 


MAXM 


The  number  of  stiffening  rings;  NR  <  3 
The  number  of  modes  with  an  initial  im- 
perfection, including  the  axisymmetric 
mode  (n  =  0) ;  INMAX  <  15 
The  maximum  number  of  Fourier  modes  to 
be  considered;  INMAX  <  MAXM  <  15 


DATA  2: 

NITMAX 

NITMIN 

MAXSTP 

KPRINT 


TFINAL 


BETA 


The  maximum  number  of  iterations  per  time 
step  (recommended  value  =  4) 
The  minimum  number  of  iterations  per  time 
step  (Recommended  value  =  3) 
The  maximum  number  of  time  steps  to  be 
computed 

Time  step  intervals  at  which  the  dis- 
placements, stress  resultants,  and  maximum 
stress  are  to  be  printed.   If  KPRINT  =  50, 
print  out  will  occur  at  the  end  of  time 
step  50,  100,  150,  etc. 
The   nondimensional  time  at  which  the 
calculation  is  to  be  terminated.   (The 
nondimensional  period  of  the  axisymmetric 
mode  is  approximately  2tt) 
Integration  parameter  of  Newmark ' s  beta- 
method  (Recommended  value  -  .1666667) 
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EPS 


Convergence  test  parameter.   Convergence 
is  assumed  when  the  magnitude  of  the 
change  in  the  acceleration  for  two 
successive  iterations  is  less  than  or 
equal  to  EPS  times  the  magnitude  of  the 
larger  acceleration.   (Recommended  value  - 
.001) 


DATA  3: 
NBC 


NSYM 


The  boundary  condition  parameter.   When  the 
shell : 

1.  Has  simple  boundary  conditions, 
NBC  =  1. 

2.  Has  clamped  boundary  conditions, 
NBC  =  2. 

Symmetry  parameter.   When  the  motion: 

1.  Is  assumed  to  be  symmetric  with 
respect  to  the  midspan,  NSYM  =  1. 

2.  Otherwise,  NSYM  =  2. 


DATA  4 


This  data  card  contains  nondimensional 
parameters  representing  the  stiffness 
and  geometric  properties  of  the  stiffen- 
ing rings.   This  card  must  be  omitted  if 
the  shell  does  not  have  stiffening  rings. 


234 


CI 


C2 


C3 


C4 


C5 


E   I  (1-v  ) 
r   x 

E   a2h  AX 


Pr  Ar 


p  hAX 


E   A  (1-v  ) 
r   r 

E    hAX 


E   I  (1-v  ) 
r   z 

E   a2h  AX 


Er  CT(l-v^) 
E  2a2h  AX 


— ,  the  eccentricity  of  the  rings 


where 


X  = 


2MM-1 


when  NSYM  =  1,  and 


DATA  5 


X  = 


MM-1 


when  NSYM  =  2. 


KR(1),  KR(2),  and  KR(3)  are  the  axial 
stations  at  which  the  rings  are  located 
(See  the  FORTRAN  listing  of  Subroutine 
WIDER1.)   If  there  are  no  rings,  omit 
DATA  5. 


DATA  6, 
DATA  7, 

and 
DATA  8: 


The  nondimensional  magnitudes  of  the 
initial  imperfection  are  placed  on  these 
cards.   For  example,  if  the  magnitude  of 
the  imperfection  in  each  mode  is  one 
one-thousandth  of  the  shell  thickness 
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WI  (1) 


WI  (2) 


WI (INMAX) 


and  the  Fourier  number  n  =  N(L),  then 

WI(L)  =  IFooxl- 

The  Fourier  coefficient  of  the  initial 

imperfection  in  the  axisymmetric  mode, 

n  =  0.   Normally,  WI(1)  =  0. 

The  Fourier  coefficient  of  the  initial 

imperfection  in  the  mode  whose  Fourier 

number  is  n  =  N(2).   (See  explanation 

of  DATA  12) 

The  Fourier  coefficient  of  the  initial 

imperfection  in  the  mode  whose  Fourier 

number  is  n  =  N (INMAX) 

If  INMAX: 

1.  Is  less  than  15,  DATA  8  must  be  omitted 

2.  Is  less  than  8,  DATA  7  must  be  omitted. 


DATA  9 : 
D 


GNU 

ALPHA 

CDAMP 


The   nondimensional  length  =  L/a.   Sub- 
routine START  sets  D  equal  to  the  non- 
dimensional  spacing  of  the  finite 
difference  stations. 
Poisson's  ratio 
-   h/a 

the  nondimensional  viscous  damping  coef- 
ficient   CDAMP  =  C  [  (1"V  ) ]  %  where  C  is 

Ep 

the  viscous  damping  coefficient  . 
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DATA  10: 
PO 
R 


(As  required  by  subroutine  PRLOAD) 
The  maximum  pressure 
The  decay  rate 
In  dimensional  units, 


-rt   „    -t/t0 


P  (t)  =  p   e    =  P   e 
o         o 


Then 


PO  =  p 


a(l-v  ) 


o    Eh 


2,  i 


2  l 


R 


S, OMEGA   -  blank.   S  and  OMEGA  are  available  to  the 
user  in  preparing  a  subprogram  to  replace 
PRLOAD  . 


DATA  12: 

N(l)  =  0 
N(2) 

N(INMAX) 


The  Fourier  numbers  of  the  initial  imper- 
fection.  These  numbers  must  be  in  one 
to  one  correspondence  with 

WI(1), ,WI(INMAX)  on  DATA  6,  DATA  7, 

and  DATA  8 .   N(l)  must  always  be  equal  to 
zero  . 


DATA  13: 
DT 


The  time  increment  to  be  used  for  the  first 
time  step  (value  used  in  calculations  of 
Chapter  V  =  .05) 


237 


DTMAX 


DATA  14: 


LAST 


The  upper  limit  on  the  time  increment. 
When  BETA  =  1/6,  DTMAX  should  not  be 
larger  than  one-fifth  of  the  shortest 
period  of  any  mode  in  the  problem. 
Normally,  the  axisymmetric  mode  has  the 
shortest  period.   (value  used  in  calcu- 
lations of  Chapter  V  =  0.7) 

This  data  card  provides  the  capability 
of  solving  several  problems  during  a 
single  computer  run. 

1.  LAST  =  ,,,,  if  one  or  more  problem 
data  sets  follow. 

2.  LAST  =  ....  if  the  present  data 
set  is  the  last  one 


Truncation  of  the  Series 

The  user  prescribes  which  modes  are  in  the  truncated 
Fourier  series  by  the  values  given  to  INMAX,  MAXM  and  N(L). 
Coefficients  having  the  same  Fourier  numbers  as 
N  (1) ,... ,N (INMAX)  are  included  in  the  truncated  series 
for  the  displacements;  and  if  MAXM  is  equal  to  INMAX,  no 
other  modes  will  be  considered  in  the  solution.   If  MAXM 
is  greater  than  INMAX,  subroutine  MODES  will  add  additional 
modes  to  the  series  for  the  displacements  until  the  total 
number  of  modes  is  equal  to  MAXM.   For  example,  if 
INMAX  =  2,  MAXM  =  4  and  N(l)  =  0,  and  N(2)  =  3,  subroutine 
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MODES  will  add  modes  six  and  nine  to  the  series  for  the 
displacements  with  zero  taken  as  the  initial  imperfection 
in  modes  six  and  nine. 

IV.   USER  ALTERATIONS  TO  OBTAIN  UNCOUPLED  SOLUTIONS 

Elimination  of  Nonlinear  Feedback  to  Axisymmetric  Mode 

In  Chapter  V  a  procedure  was  outlined  for  identifying 
the  Mathieu  modes.   The  procedure  involves  computing  the 
response  of  the  flexural  modes  whose  frequencies  are  close 
to  one-half  the  frequency  of  the  axisymmetric  mode.   In 
setting  up  the  program  for  such  an  identifying  run,  the 
nonlinear  terms  in  the  equations  of  the  axisymmetric  mode 
are  set  equal  to  zero.   The  following  alterations  of  the 
program  will  eliminate  the  appropriate  nonlinear  terms. 

In  subroutine  LAM1N2 : 
After  the  statement 

11    D(j)  12   L  =  1,  MNMAX0 
insert 

G(J)  TCj)  2  0 
10  0  CONTINUE 

The  second  inserted  statement  is  necessary  to  avoid  an 

error  message  from  the  Fortran  compiler. 

In  subroutine  TNETA: 
After  the  statement 

9  D(J>  10    L  =  l,MNMAX(j) 
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insert 

G§    T(|)  11 
100  CONTINUE 
The  second  inserted  statement  is  necessary  to  avoid  an 
error  message  from  the  Fortran  compiler. 

Elimination  of  Nonlinear  Coupling  Between  the  Flexural  Modes 

In  Chapter  V  a  comparison  was  made  between  a  fully 
coupled  nonlinear  solution  and  one  with  limited  coupling. 
To  eliminate  the  intermodal  coupling  between  the  flexural 
modes,  but  yet  retain  the  coupling  between  each  flexural 
mode  and  the  axisymmetric  mode,  the  following  alterations 
are  required: 

In  subroutine  LAM1N2 : 
After  the  statements 

IF(N(M)  .EQ.O)  G(|)  T(J)  11 

MAXL  =  MAXS(M) 
insert 

MAXL  =  1 
After  the  statement 

8  MAXL  =  MAXD(M) 
insert 

MAXL  =  1 
After  the  statement 

10  IF  (MAXSY(M)  .EQ.O)  G(|)  T<\)    13 

insert 

G(j)  T(j)  13 
100  CONTINUE 
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The  second  inserted  statement  is  necessary  to  avoid  an 
error  message  from  the  Fortran  compiler. 

In  subroutine  TNETA: 
After  the  statements 

IF  (N(M)  .EQ.O)  G(j)  T(J)  9 

MAXL  =  MAXS (M) 
insert 

MAXL  =  1 
After  the  statement 

6  MAXL  =  MAXD (M) 
insert 

MAXL  =  1 
After  the  statement 

8  IF(MAXSY(M) .EQ.O)  G0  T0  11 

insert 

G0  T<|>  11 
100  CONTINUE 

The  second  inserted  statement  is  necessary  to  avoid  an 

error  message  from  the  Fortran  compiler. 

V.   STORAGE  REQUIREMENTS  AND  EXECUTION  TIME 

The  program  requires  approximately  125K  bytes  of 
storage  space  in  the  IBM  360  Model  67  memory  when  sub- 
routines AMPMAX,  WI1SIN,  and  WI2SIN,  or  WIDERl  and  WIDER2 
are  removed  from  the  deck.   This  figure  includes  the  space 
required  for  the  execution  of  the  required  input  and  out- 
put routines. 
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The  approximate  execution  time  for  calculating  a 
completely  coupled  nonlinear  solution  until  t  equal  two- 
hundred  is  shown  in  Figure  D.2  as  a  function  of  the  number 
of  modes  in  the  Fourier  series. 

VI.   PROGRAM  LISTING 

The  listing  of  the  program  is  contained  on 
pages  244  -  279. 
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